The occurrence of seismicity accompanies mining activity conducted primarily in underground mines that extract various kinds of raw mineral material. This is true for numerous mining regions worldwide, where seismicity often becomes the cause of various geodynamic phenomena, particularly rock bursts as well as vibrations of the Earth’s surface (Konicek et al. 2019; Feng 2017; Li et al. 2007; Mendecki 2015).

It should be stressed that rock bursts result in the loss of working functionality to a diverse degree, as well as in machine and equipment damage, and accidents or disasters among the personnel, thereby constituting a major technical, social and economic problem (Drzewiecki and Piernikarczyk 2017; Li et al. 2007; Brady and Brown 2007). Furthermore, major mining-induced tremors lead to various repercussions on the surface in the form of residential building and technical infrastructure damages (Dubiński et al. 2020; Kalab 2018).

An example of such a situation is the Upper Silesian Coal Basin (USCB), where extraction has been conducted for over 200 years by both Polish and Czech underground hard coal mines (Dubiński et al 2019; Konicek et al. 2019). The evolution of seismicity and rock bursts in Polish USCB mines over the period of 1980–2020 is presented in Table 1. The data included in the table is limited only to stronger phenomena with a seismic energy E_S ≥ 10⁵ J, which corresponds to a local magnitude of about M_L = 1.6.

Despite the significant decrease in coal extraction, a proportional decrease in seismicity has not been observed. Currently, seismic hazards can be found in nearly all the mines, and can thus be considered a common hazard. For this reason, the issue of mining-induced seismicity and rock bursts has always been one of the priorities for research conducted as part of various types of projects. Furthermore, all mines that conduct mining activity in coal seams at risk of rock bursts are equipped with modern seismological, seismoacoustic and seismic apparatus that forms a geophysical rock mass monitoring and measuring system. A daily seismic hazard evaluation is established for each mining area exposed to the hazard, based on the results of constant observation and direct underground measurements, which forms the basis for rational decision-making within the scope of occupational safety and the required preventive measures (Dubiński and Konopko 2000).

Such a state necessitates the constant enhancement and development of new measurement methodologies with the purpose of improving the quality of the daily seismic and rock burst hazard evaluations.

The presented attempt at correlating selected seismological data and underground seismic measurement results indicates the possibility of increasing the seismic hazard evaluation precision. The obtained results and the proposed new comprehensive evaluation are based on tests conducted in hard coal mines in the Upper Silesian Coal Basin, which are characterised by the occurrence of major seismic hazards.

The testing ground encompassed six longwalls under extraction by two mining plants belonging to Polish Mining Group Ltd. (Polska Grypa Górnicza S.A.). There were four longwalls located in the coal seam 504 (longwalls 002 and 003) and 506 (longwalls 1 and 2) of the Bielszowice coal mine. The second mine was Myslowice-Wesola mine where two longwalls were located in the coal seam 510 (longwall 01Aw and 02Aw).

The both mines belong to mining plants where seismic and rock burst hazards have constituted the most significant natural hazards for many years. The characteristic geological and mining conditions occur in both mines, which include (Dubiński et al. 2019):

1. (1)

   Great mining depth, often exceeding even 1000 m, resulting in high geostatic stress,

2. (2)

   The presence of thick, high-strength rock strata in the roof rock,

3. (3)

   Coal seams propensity for rock bursts,

4. (4)

   Complex tectonic structure in the form of faults and other geological disturbances,

5. (5)

   High extraction intensity resulting in the considerable vertical range of the zone of rock mass equilibrium disturbance,

6. (6)

   The multiple-seam character of the deposit and the many years of its extraction, resulting in major rock mass damage and the presence of coal seam and face remnants in other coal seams, leading to the formation of stress concentration zones.

Both the mines are equipped with digital seismological systems with software and computer bases as well as seismic apparatus for conducting underground seismic measurements.

It should be stressed that the mines where the testing grounds were located are found in different parts of the USCB, therefore the obtained test material was diverse. This was a necessary condition for testing and calibrating the measurement parameters that constitute the basis of the new comprehensive method.

The comprehensive character of the developed seismic hazard evaluation method consists in the application of two different components:

1. (1)

   Mine seismology utilising the results obtained by the mine measurement network,

2. (2)

   Seismic method where the measurement range encompasses the planned longwall panel.

Each of these component methods makes it possible to determine the parameters characterising the rock mass from the perspective of its propensity for seismicity. The first parameter was b coefficient of the Gutenberg–Richter distribution (Aki 1965) and the second was the maximum value of longitudinal seismic wave \({V}_{P\mathrm{max}}^{\mathrm{roof}}\) propagating velocity in roof strata.

3.1 Seismological method—the Gutenberg–Richter law

In the case of mining-induced seismicity, similar as to earthquake-related seismicity, it is described by the relationship between the number of the generated tremors and their scale expressed in the magnitude or the seismic energy parameter. This relationship is known in literature as the Gutenberg–Richter law (GR), in the following form (Aki 1965; Utsu 1965):

It is also used in mining seismology (Xu et al. 2014; Mutke et al. 2016).

For mining-induced tremors, instead of M_L the applied parameter is the tremor seismic energy E_S, which for the seismicity occurring in the USCB is related to the local magnitude by the following relationship (Dubiński and Wierzchowska 1973):

where N is the number of tremors in a magnitude class (M, M + dM), a and b are constant coefficients, M_L is local magnitude, E_S is seismic energy, J.

The a and b coefficients in Eq. (1) have their own physical interpretation. Particularly the b coefficient provides valuable information regarding the rock medium where the stress and strain processes take place and the foci of the generated tremors are formed. The reason for this is that this coefficient describes the relationship between strong and weak seismic phenomena. A high b value indicates the presence of a greater number of energetically weak seismic phenomena in a given tremor group, whereas low b values correspond to an opposite relationship, i.e. to the presence of a greater number of tremors with higher seismic energy (Gibowicz and Kijko 1994; Mutke et al. 2016).

Therefore the b coefficient has great significance in the development of the new seismic hazard evaluation method. That is because its decreasing value indicates an increase in stress within the focus, whereas the stress falls as the value increases. These variations, expressed as the anomaly of the b value registered during longwall extraction, signify that the seismic hazard intensifies or decreases and make it possible to evaluate its level.

3.2 The seismic geotomography method

The method based on seismic measurements conducted in underground workings is commonly applied in Polish hard coal mines. The object of the measurements is the structure of the coal seam and the surrounding rock, which undergoes changes over time as a result of the developing stress and strain processes generated by the conducted extraction (Maxwell and Young 2012; Chen et al. 2015). These variations are reflected in the seismic wave propagation velocity in the aforementioned rock media, most frequently the longitudinal wave, both with regard to the absolute value of the velocity parameter and the occurring its anomalies (Dubiński and Konopko 2000).

Years of experience in utilising the seismic method in solving various geomechanical and geological problems found in hard coal mines demonstrate its particular usefulness in the pre-emptive identification of seismic hazard zones. In practice, the method fulfils a role similar to the method for recognition the geological and mining conditions determining the level of the seismic and rock burst hazards. It should be stressed that there are many variants of the measurement methodology, selected and applied depending on the local geological and mining situation and the specifics of the problem to be solved (Cai et al. 2015; Dou et al. 2012).

As the second component of the comprehensive method, the seismic method is based on the results of seismic longwall panel ray scanning performed periodically both before extraction is commenced and during its conduction. The seismic wave receivers are geophone probes installed in holes drilled in the side of the gateroads.

Meanwhile the seismic wave excitation process utilises minor explosive charges. The distribution of both the excitation and the seismic wave reception points is planned in such a way so as to obtain the best possible coverage of the coal seam part scanned by the seismic rays. The interpretive process utilises geotomographic wave field reconstruction procedures, most often referenced to the longitudinal seismic wave propagating in the roof rock (Dubiński and Konopko 2000; Nolet 2008; Young and Maxwell 1992). The result of this reconstruction is a velocity isoline distribution map for the particular seismic wave in the scanned longwall panel area. It serves as the basis for calculating the value of the seismic parameter \({V}_{P\mathrm{max}}^{\mathrm{roof}}\) which is the second component of the comprehensive seismic hazard evaluation method.

The developed methodology encompasses a number of basic procedure steps concerning both the b value of the GR law, and the second parameter of the maximum velocity \({V}_{P\mathrm{max}}^{\mathrm{roof}}\) determined based on the velocity isoline distribution map.

1. (1)

   Determining the threshold value of the magnitude \({M}_{\text{L}}^{t}\) or seismic energy \({E}_{\text{S}}^{t}\) based on sets of tremors registered in a given extraction area by analysing the number of tremors in individual energy intervals. The threshold value signifies the lowest value of the magnitude or seismic energy above which the GR relationship is true (the determination scheme is displayed in Fig. 1.

2. (2)

   Selecting the time scale for calculating the b value (the time scale should include at least 20–30 tremors).

3. (3)

   Determining the calculation step (as the seismic hazard evaluation is performed for each successive day, the adopted step is equal to 1 day).

4. (4)

   Calculating the b value using the maximum likelihood method and a standard deviation σ_b based on the relationships provided by Aki (1965) and Utsu (1965):

   $$b = \frac{\log e}{{\overline{M} - M_{\min } }}$$

   (3)
   $$\sigma_{b} = 2.3b^{2} \sigma_{{\overline{M}}}$$

   (4)

   where \(\sigma_{{\overline{M}}} = \sqrt {\sum\nolimits_{i = 1}^{n} {\frac{{\left( {M_{i} - M} \right)^{2} }}{{n\left( {n - 1} \right)}}} }\), e is the Euler’s number; \(\overline{M }\) is average tremor magnitude in a set, M_min is minimum threshold magnitude, M_i is magnitude in the i-th calculation.

5. (5)

   Determining the average b_med value, calculated from the date of longwall extraction commencement to the date of calculation performance.

   $${b}_{\mathrm{med}}=\frac{\left({b}_{1}+\dots +{b}_{n}\right)}{n}$$

   (5)

   where b₁, …, b_n—b values during individual days of extraction.

6. (6)

   Calculating the Gutenberg–Richter anomaly zA_GR using Eq. (6).

   $$z{{{A}}}_{{{G}}{{R}}}=\left[\frac{\left({{{b}}}_{\mathrm{med}}- {{b}}\right)}{{{{b}}}_{\mathrm{med}}}\right]\bullet 100\boldsymbol{\%}$$

   (6)

   where b_med is average b value calculated since the date of longwall extraction commencement.

Example magnitude threshold value determination

The computer program GMB was developed for the practical application of the presented methodology and the daily zA_GR parameter calculations, and it currently finds common use in mining plants.

1. (7)

   Assigning the appropriate weights to the zA_GR parameter values, quantifying the seismic hazard level forming in the zone of the developing seismic activity in the rock mass, as displayed in Table 2.

Table 2 presents two situations differing in the relationship between the current b value and its average b_med value as well as the criterion value b = 1.5 adopted as the threshold above which seismic hazards are not present. When such relationships occur, it is not necessary to analyse the GR anomaly zA_GR. However, when b < 1.5, a seismic hazard can occur only when the condition b < b_med is fulfilled at the same time. In such cases the zA_GR parameter is analysed, and three hazard levels are distinguished based on its value, with assigned weights of 1, 2and 3.

1. (8)

   Determining the maximum seismic wave propagation velocities \({V}_{P\mathrm{max}}^{\mathrm{roof}}\) based on an isoline map of longitudinal wave propagation velocity in the roof rock, obtained as a result of seismic ray scanning. These are determined in specific fragments of the longwall panel before the longwall face, whose locations vary based on the rate of longwall advance.

2. (9)

   Assigning the appropriate weights to the individual \({V}_{P\mathrm{max}}^{\mathrm{roof}}\) parameter values, which quantify the roof rock propensity for tremor generation, based on the application of an existing seismic scale (Dubiński and Konopko 2000), presented in Table 3.

The values of weights for the b coefficient of the Gutenberg–Richter distribution and the longitudinal wave velocity in the surrounding rocks \({V}_{P\mathrm{max}}^{\mathrm{roof}}\) were selected empirically based on the experience related to the actual seismic hazard observed for each longwall. The weights of both parameters were determined in such way that the obtained value for each of them and then the final value being the sum of zA_GR and z\({V}_{P\mathrm{max}}^{\mathrm{roof}}\) were adequate to the seismic hazard observed during the longwall operation. Various weights were checked many times to reflect the seismic hazard as accurately as possible.

(10) Determining the resultant seismic hazard level Z_GR+SG according to the comprehensive method based on the G-R law and the results of seismic geotomography SG. The values of the numerical intervals of the sum of the weights used to prepare the seismic comprehensive hazard assessment Z_GR+SG were also selected empirically based on the observations made for each of the longwalls.

The sum of weights obtained from the constituent methods classifies an area of conducted longwall extraction as exposed to one of four different seismic hazard levels. Table 4 presents the final evaluation scale of the seismic hazard Z_GR+SG based on the two aforementioned component methods.

It should be emphasized that all weight values for the parameters zA_GR and z\({V}_{P\mathrm{max}}^{\mathrm{roof}}\) as well as the seismic hazard assessment Z_GR+SG were determined for the set of six analyzed longwalls operated in the USCB, in coal seams 504, 507 and 510, at depths from about 800 m to about 1050 m.

5.1 Geological and mining conditions

Longwall 02Aw was another longwall excavated in the first top slice, of 3 m high, of the coal seam 510 in the A section, the thickness of which is variable and in this region varies from 7.7 to 10.7 m. The roof of coal seam 510 exhibits sandstone and shale strata with a uniaxial compressive strength of about 60 MPa. These rock strata dip towards the S-W at an angle of about 6°. The mining depth varies from 800 to 850 m. Coal seam 510 in the longwall area was classified as exposed to a second level rock burst hazard. Figure 2 presents a fragment of the coal seam 510 map, with the location of the longwall 02Aw panel and the geological profile of the roof strata in the area.

Map of coal seam 510 with the location of longwall 02Aw and the geological profile of the roof strata

The extraction area was characterised by a complex faulting tectonic structure, and coal seam 510 was mined in the hanging walls of two major faults with throws of 175 m and 250 m. In front of the longwall 02Aw occurred the edges of the above-lying coal seams with numbers 401, 404/1, 404/5 and 405/2, deposited at distances from 175 to 260 m.

The mining plant applied active preventive measures in the mined plot of coal seam 510, in the form of torpedo and destressing blasting for rock burst hazard prevention. The area of the studied longwall 02Aw was under constant seismic monitoring by a well-advanced mine seismological network, which made it possible to locate tremor foci with an average accuracy of 20 m.

5.2 Characteristics of seismicity during the mining of longwall 02Aw

Longwall 02Aw was the second successive longwall driven in part Aw of coal seam 510. The seismic activity developed together with the longwall face advance and the increase in the goaf area. Table 5 presents the energy and quantity characteristics of the registered seismicity associated with this longwall.

This seismicity should be classified as medium, where the number of high-energy tremors (E_S ≥ 10⁵ J) is low and their occurrence is sporadic. The focus epicentre map for tremors with seismic energy E ≥ 10³ J that occurred during the extraction of longwall 02Aw is presented in Fig. 3.

Distribution of seismic epicentres in the longwall 02Aw panel

It should be noted that the registered seismicity was nearly twice as great during the extraction of longwall 02Aw compared to the previously mined longwall 01Aw. 904 tremors were registered in total. A tremor with a high seismic energy E_S = 3 × 10⁸ J occurred during the final stage of the extraction, and was accompanied by a rock burst.

5.3 Underground seismic measurements and obtained results

The first cycle of seismic ray scanning was performed on a 300-m-long section of the longwall 02Aw panel before the extraction of this longwall was commenced. The second cycle, however, was carried out after 5 months, when the longwall face was 440 m before its end. Figures 4a and b present the longitudinal seismic wave propagation velocity isoline distribution maps obtained as a result of these two measurement cycles.

Results of seismic measurements in the longwall 02Aw panel using active seismic geotomography: a 1st cycle b 2nd cycle

Based on the results of the first seismic measurement cycle, it was concluded that the longitudinal seismic wave propagation velocity in the surrounding rock varied within 3480 m/s to 4260 m/s (average velocity 3950 m/s). Therefore, the maximum velocity values for seismic waves propagating in the roof rock of the first part of the longwall encompassed by the measurements have been classified as having a medium propensity for generating seismicity. During the second measurement cycle, the longitudinal seismic wave propagation velocity variations in the roof rock ranged within 3470 m/s to 4150 m/s (average velocity 3820 m/s). Therefore the roof rock classification for the second part of the longwall panel was similar, with the exception that lower longitudinal wave velocities were found over a significantly greater scanning area, which indicates a lower propensity for generating seismicity.

Local seismic anomalies relative to the above average velocities can be located in both the velocity isoline maps. They indicate the presence of zones of increased stress as well as destressed zones, generated as a result of prior preventive measures taken against rock bursts.

5.4 Gutenberg–Richter (GR) law parameter determination results

Medium seismicity was found during the extraction of the first, 300 m long part of longwall 02Aw (part I), as 298 tremors were recorded, with seismic energy ranging from 10² to 7 × 10⁴ J. The determined threshold of a minimum magnitude \({M}_{\text{L}}^{t}\) = 0.63 and corresponds to seismic energy threshold \({E}_{\text{S}}^{t}\)=10³ J. The b value calculations were carried out according to the methodology provided, in 20-day time scales with steps of 1 day. The calculations were performed using the GMB program. Variations of the b value and its average value b_med together with a histogram of the maximum seismic energy values during individual days of extraction are presented in Fig. 5a.

Variations of the b and b_med parameters and \({E}_{\text{S}}^{\text{max}}\) along the longwall 02Aw panel in coal seam 510: a part I b part II

On the other hand, the further extraction of longwall 02Aw in the second, 440-m-long part of the longwall panel (part II), where the longwall face approached the final section by the boundary of goafs belonging to the prior extracted part of coal seam 510, was characterised by significantly higher seismicity. 425 tremors were registered, including phenomena with seismic energy values higher than 10⁵ J. The calculated b value for the analysed period underwent a significant decrease. Figure 5b presents the variations of the b value and its average value b_med together with a histogram of the maximum seismic energy values \({E}_{\text{S}}^{\text{max}}\) recorded during successive days of extraction.

5.5 Correlation of the seismological and seismic results

The correlation between the isoline distribution of the longitudinal seismic wave propagation velocity in the roof rock and the registered seismicity described by the b value during the extraction period in part I confirms that low seismicity corresponds to zones with lower values of longitudinal seismic wave propagation velocity in the roof rock, whereas the b and b_med parameters in such cases exhibit growing trends. On the other hand, when the longwall 02Aw face started approaching zones with increased longitudinal wave velocities, the b value began to decrease rapidly, though the simultaneous reduction of the b_med parameter was much slower.

However, during the extraction in part II, the area of the longwall 02Aw goafs began to increase, which resulted in the occurrence of tremors with greater seismic energy values in the seismicity image. This led to a significant reduction of the b value, below the criterion value of 1.5 and also below the value of b_med. Such a relationship between these values corresponds to a higher seismic hazard level. At the same time, a clear difference can be observed in the longitudinal wave velocity isoline distribution associated with the roof rock. There are significantly lower longitudinal seismic wave velocities in the roof rock along the majority of the scanned panel, which indicates that the roof rock propensity for generating seismicity is rather low.

An advantage of the comprehensive seismic hazard evaluation using the method in question is that its constituents exhibit different reactions to certain factors that determine seismicity. For this reason, the obtained results are closer to the actual state. Such a situation can be found in the extraction of longwall 02Aw panel part II, where the seismicity increases (as indicated by the seismological component and the b value), but the second component alleviates this increase, thereby revealing it to be at a medium level rather than at a very high one. The tremor with an energy of 3 × 10⁸ J was a phenomenon of a regional character that occurred in an area of significant overlying roof rock stratum strain as a result of prior extraction and the complex faulting tectonic structure.

5.6 Comprehensive seismic hazard evaluation

Analysing the results obtained by the individual constituent methods (b value and velocity \({V}_{P\mathrm{max}}^{\mathrm{roof}}\)) indicates that the components exhibit various contribution to the resultant seismic hazard evaluation. The longitudinal seismic wave propagation velocity in the roof rock ranges within 3500–4500 m/s, which entails a medium propensity for generating seismic phenomena in the entire mined longwall 02Aw panel in coal seam 510. Therefore, the G-R law's b value determines the seismic hazard level in this case both in terms of its increase and decrease.

Areas where the basic condition confirming the presence of seismic hazard zones is fulfilled (b < b_med and b < 1.5) can be discerned in Fig. 5a and b. A significant difference in the size of these zones between part I and part II should be noted. In part I, the zone is small and encompasses a time of about 10 days of extraction, whereas in part II the above condition for seismic hazard presence is fulfilled during the entire extraction period of the part in question. In the first narrow hazard zone in part I, the comprehensive evaluation indicated a low hazard (level b), with no hazard (level a) for the remainder of the extraction period. On the other hand, the significant majority of part II was exhibited a weak hazard (level b), but a medium hazard (level c) occurred several times, which indicated a rising seismic hazard. Among other factors, it was such an increase that transpired before the occurrence of the high-energy tremor with E_S = 3 × 10⁸ J that was accompanied by a rock burst.

Per the assumptions of this method, it is intended to provide an auxiliary measure of seismic hazard evaluation together with the standard set of other geophysical methods that have found application in this field thus far.

5.7 Comprehensive method validation

When testing and calibrating new measurement methods for their future application in practice, carrying out their validation, in this case with regard to increasing the quality of the seismic hazard evaluation, is an essential operation. The reason for this is that the effect of implementing such a method can be evaluated based on the validation. Under the conditions of a mine geophysical station that conducts constant monitoring with the purpose of providing a daily seismic hazard level evaluation based on applicable instructions within this scope, the best point of reference is to assess the conformity of the daily seismic hazard level evaluations. Seismic hazard level quantification is conducted based mainly on the parameters: maximum value of seismic energy and a value of cumulated energy related to 5 m longwall advance. Additionally such parameters as energy index EI and PPV (Particle Part Velocity) may be used (Mutke et al. 2015).

For this purpose, the evaluation results obtained by the two methods have been compiled in Table 6 for a set of tremors with seismic energy E_S ≥ 10⁴ J, recorded during the extraction of longwall 02Aw in coal seam 510.

A comparative analysis of the seismic hazard evaluation results obtained by two different methods—the reference method and the new tested method—indicates the following conformity (Table 7).

There were 24 conforming evaluations during a given days, which constitutes 42% of the total number, as non-conforming evaluations amounted to 33 (58%). A discrepancy in the medium seismic hazard level evaluation can be observed, where the evaluation per the new method exhibits a significantly greater number of indications. On the other hand, the reference method exhibits a greater number of low level evaluations.

The overall assessment of the new method in the light of these results is positive. Conformity on a level of 42% is not an unsatisfactory result in this case, given that the two methods—the new and the reference method—are based on parameters that exhibit different reactions to the development of strain and stress processes determining the propensity for generating tremors by the rock mass. The new method, utilising the b value, is physically more sensitive to the variations in the rock mass structure as a result of the aforementioned processes. This may be the reason for the more frequent indication of higher seismic hazard levels. Furthermore, full conformity between the indications of the two methods would not provide any new information regarding the complex processes determining the seismicity of the rock mass and the risks it poses to the activity of mining plants.

Nevertheless, the conformity assessment results of the methods in question indicate the need to:

1. (1)

   Verify the criterion values for the individual hazard levels in the methodology utilising the GR law's b value, which ought to be calibrated with reference to each new area where the method is applied,

2. (2)

   Consider a broader use of the isoline distribution maps for the longitudinal seismic wave propagation velocity in the roof rock (using only the maximum velocity value may not fully reflect the variability of this parameter in the scanned longwall panel).

The assumption that, at its current stage of development, the new comprehensive seismic hazard evaluation method based on the GR law's b value and the longitudinal seismic wave velocity parameter obtained using seismic geotomography will become an auxiliary method for the other currently applicable methods is fully justified. Its improvement will be enabled by further testing and evaluation parameter calibration.

1. (1)

   The common application of geophysical measuring systems with digital registration in Polish hard coal mines introduced new possibilities for data processing as well as for using new parameters that characterise the variations occurring in the rock medium structure under the influence of the conducted coal seam extraction.

2. (2)

   The conducted tests confirmed that the comprehensive solution using the b parameter of the Gutenberg–Richter law and the longitudinal wave velocity distribution in roof strata provides such information.

3. (3)

   The proposed new mine seismic hazard evaluation criteria can expand the practical capabilities of the commonly applied measurement methods.

4. (4)

   The seismic hazard evaluation credibility validation for the new comprehensive evaluation method indicates the presence of a discrepancy relative to the currently applicable method based on the seismic energy of the registered tremors.

5. (5)

   It is recommended to continue the research on improving the developed solution, particularly with the purpose of verifying the applied criterion values for the seismic hazard evaluation.