Dongtan Coal Mine is located in Yanzhou, Shandong, China, and attracted widespread attention due to frequent strong mine tremors (Figure 1A). In the course of mining operations in the Dongtan No. 6 mining area, a multitude of tremors were detectable by the earthquake monitoring network. These tremors resulted in the shaking of residential structures in proximity to the mine, thereby engendering significant social repercussions. The current main coal seam is 3^upper with a buried depth of about 680 m. The thickness of 3^upper coal seam is in the range of 4.12–6.79 m, with an average value of 5.41 m. The coal seam is nearly horizontal and the geological structure is relatively simple. The maximum horizontal principal stress is 24.96∼27.12 MPa at the direction of SE 148.93°∼150.00°. The maximum horizontal principal stress is 2.15∼3.29 times the minimum horizontal principal stress and 1.44∼1.72 times the vertical stress.

Dongtan coal mine is divided into 7 mining areas, and the No.6 mining area is located in the south of the coal mine. Longwall (LW) fully mechanized mining technology is adopted. Four working faces have been exploited with the mining sequence of LW63^upper04 ⇒ LW63^upper05 ⇒ LW63^upper03 ⇒ LW63^upper06 (Figure 1B). The former three are finished and turned into goaves. As of October 2022, the LW63^upper06 has been advanced to 879 m (the panel length is 1,500 m). The lithology of the roof and floor of the 3^upper coal are mainly sandstone with high strength, and their main mechanical parameters have been measured in the laboratory (see Supplementary Table S1).

A 16-channel microseismic (MS) monitoring system, developed by the Institute of Innovative Technologies EMAG, Poland, has been installed to determine the location and energy of mine tremors. DLM microseismic sensors with an operating frequency of 0.1–200 Hz are included in this system. MS signals are transmitted via electric transmission lines to the surface. Eight sensors (4#, 5#, 6#, 14#, 15#, 17#, 28#, 29#) are installed in different positions of underground roadways in the No. 6 mining area and two sensors (16#, 29#) are installed on the ground to form a joint stereoscopic MS monitoring network (Figure 2A), which improves the positioning accuracy of the microseismic monitoring system in the vertical direction. The sampling frequency range of these sensors is 500 Hz. The No. 6 mining area of Dongtan Coal Mine experienced a substantial number of strong mine tremors during previous coal mining. These tremors prompted considerable social anxiety and drew significant government attention. Specifically, 15 mine tremors, with magnitudes exceeding 1.5, were recorded on LW63^upper03, while LW63^upper04 and LW63^upper05 reported 39 and 56 such tremors, respectively. Compared with other working faces, mine tremor on working face LW63^upper06 is more severe. Since the start of mining in February 2020, there have been more than 34 large energy mine tremors with a magnitude of more than 2.0 (Figure 2C). It can be seen from Figure 2B that strong mine tremors mainly occur in the roof rock layer, especially near the lower boundary of the “red layer”.

In Dongtan Coal Mine, there have been numerous strong mine tremors, but they rarely lead to roadway rock bursts. This phenomenon might be closely linked to the vibration characteristics of these tremors and the distance between their source and the roadways. Figure 3 presents waveform diagrams of two distinct types of strong mine tremors. Wave 1, originating from the red layer located 180 m away from the coal seam, is classified as a far-field wave. Wave 2, originates from the central sandstone layer just 30 m from the coal seam. Both waves are recorded by the 14# microseismic sensor. Both mine tremors display approximately 2 s of intense fluctuations, followed by a phase of minor amplitude variations. This temporal pattern sets them apart from natural earthquakes, which typically last for over 10 s. Notably, Wave 2 has a higher dominant frequency than Wave 1. Additionally, the maximum amplitude of Wave 2 seems to exceed the detection range of the microseismic sensor. Zhou et al. (Zhou et al., 2016) used numerical simulations to show that medium properties cause frequency attenuation in blast waves over distance. Similarly, the stark contrast between the two waveforms might result from either frequency and amplitude attenuation during propagation or significant differences in tremor sources, making a clear explanation challenging. It is worth noting that what we analyze here is the frequency and vibration speed when the vibration propagates from the source to the roadway.

Ambient noise and instrument-related interference can complicate the recorded waveforms. Variational Mode Decomposition (VMD), a non-recursive signal processing method developed in 2014 was utilized to address this (Dragomiretskiy and Zosso, 2013; Liu et al., 2016; Liu et al., 2023). VMD is renowned for its effective noise reduction capabilities and its ability to mitigate the issue of frequency aliasing commonly seen in empirical mode decomposition (EMD). Figure 4 shows the result of two typical waveform’s VMD and Hilbert spectrum. The main frequency of wave 1 is in the range of 0–5 Hz, and the duration is longer than 2 s (Figure 4B). The components of wave 2 are quite complex, with frequencies ranging from 0 to 80 Hz, but the components in the range of 20–80 Hz last for about 1 s and disappear, and the components below 10 Hz last for more than 3 s (Figure 4D). The analysis above indicates that the high-frequency component of the waveform decays much faster than the low-frequency component.

To analyze the effect of vibration on the roadway after it is generated from the tremor source, a numerical model of fine calibration is established for analysis. According to the actual measured borehole histogram in study area, a global model was established after proper simplification of the rock strata (Figure 5A). The model size was 150 m ×30 m × 164 m (length × width × height). The number of zones in numerical modeling is 355,000, the size of the rock strata grid is less than 2 m, and the size of the coal seam grid is about 0.5 m. The four sides and the bottom of the model were constrained by velocity. The average burial depth of LW63^upper06 is 680 m, so vertical stress of 17.2 MPa was applied to the top of the model. Horizontal stress is applied by the real ground stress test results detailed in section 2.1. The Strain-softening model was used to simulate the coal seam, and the Mohr–Coulomb model was used to simulate other rock strata.

The simulation process, detailed in Figure 5B, began with the assignment of initial parameters and equilibrium calculations. Subsequently, the roadway was excavated, and the model was recalculated to reach equilibrium. Finally, the dynamic calculation model was opened, various dynamic loads were applied, and the study examined the propagation and attenuation patterns of these dynamic loads and their rockburst potential on the roadway.

The physical and mechanical parameters of the rock mass were determined through laboratory tests conducted on test specimens, combined with a parameter conversion method previously proposed. To establish this relationship, data from 44 different sources were collected to compare experimental mechanical parameters with simulated ones, resulting in a linear conversion relationship as follows (Mohammad et al., 1997): y μ = 1.0556 × x μ (1) y E = 0.4692 × x E (2) y U T S = 0.4943 × x U T S (3) y U C S = 0.2837 × x U C S (4)

Where, y μ and x μ represent simulation and laboratory Poisson’s ratios, respectively; y E and x E represent simulation and laboratory Young’s modulus, respectively; y U T S and x U T S represent simulation and laboratory Uniaxial tensile strength, respectively; y U C S and x U C S represent simulation and laboratory Uniaxial compressive strength, respectively.

Dongtan Coal Mine extensively measured the mechanical parameters of coal and rock (see Supplementary Table S1). Utilizing the conversion relationships outlined in formulas (1), (2), and (3), the mechanical parameters for the model’s rock strata were determined (shown in Supplementary Table S2). The coal seam employs a strain-softening constitutive model, with its mechanical parameters calibrated through numerous iterations.

As depicted in Figure 7, a total of 12 monitoring points were strategically positioned within the meticulously calibrated numerical model to meticulously capture and assess alterations in vibration velocity induced by dynamic loading. This arrangement comprises nine monitoring points distributed along the path from the source to the roadway’s apex, with an additional three points situated at the roadway’s base. In alignment with field conditions, the dynamic load source is applied to the primary key layer, known as the “red layer.” The figure also delineates the P-wave’s propagation pattern, encompassing a 10 m×10 m spatial extent, with its initial vibration directed forward along the z-axis.

In the Mohr-Coulomb constitutive, the plastic zone represents material failure, providing insight into the extent of material damage. This section delves into the influence of dynamic loads on the roadway by examining the plastic zone under different dynamic loading conditions.

Brittle shear ratio (BSR) is usually used to evaluate the rockburst potential in coal and rock mass, which is expressed as (Vennes and Mitri, 2017): B S R = σ 1 − σ 3 U C S (6) where σ 1 is maximum principle stress, σ 3 is minimum principle stress, and UCS is uniaxial compressive strength.

Based on the measured UCS of coal and rock (see Supplementary Table S1), combined with the σ 1 and σ 3 obtained through simulation, the BSR value can be determined. The BSR criteria for estimating the coal and rock mass damage and rockburst potential are shown in Table 1 (Castro et al., 2012). It can be estimated that the rock burst potential is major if the BSR value is higher than 0.7.

Figure 19 provides a comprehensive synthesis of previous section findings, allowing for a comparative analysis of P-wave and S-wave propagation characteristics. As shown in Figure 19A, the attenuation coefficient of P-wave initially decreases, with subsequent minor fluctuations, as the dynamic load frequency increases. The maximum attenuation coefficient values are observed at dynamic load frequencies of 2 Hz and 5 Hz. Conversely, S-wave attenuation coefficients show an initial increase followed by a decrease, stabilizing at higher frequencies, with peak values at 10 Hz and 20 Hz. Notably, S-wave attenuation coefficients are consistently lower than those of P-wave. This conclusion is consistent with that obtained by Torres using seismic observation data in the northern Gulf of California and Mexico (Bischoff et al., 2010).

As shown in Figure 19B, with the increase of frequency the peak vibration velocity induced by P-wave experiences initial minor fluctuations, followed by a rapid increase, and ultimately stabilizes at a high level. Conversely, the pattern for peak vibration velocity caused by S-wave is simpler, showing a consistent increase with rising frequency. However, there is a distinctive spike in peak vibration velocity at 30 Hz. This may be due to the resonance of the roadway caused by the S-wave at 30 Hz, resulting in severe plastic failure of the roadway. This is consistent with the results of Liu et al., who discovered that the resonant frequency of the roadway in Zhuji Coal Mine is 35 Hz (Liu et al., 2019). The trend for the amplification coefficient, both for P-wave and S-wave, closely follows that of the peak vibration velocity. Notably, beyond a frequency of 50 Hz, the amplification coefficient for S-wave significantly surpasses that of P-wave, highlighting a growing disparity between them with increasing frequency. This is consistent with the law of the plastic zone volume in Figure 21 below, where the total plastic zone volume caused by the S-wave is higher than that of the P-wave.

The 30 Hz S-wave exhibits distinct characteristics that merit further analysis, as illustrated in Figure 19B. Figure 20A clearly shows that the 30 Hz S-wave induces a higher vibration velocity at the P9 compared to other frequencies. Additionally, we have delineated the apparent fluctuation segment of each velocity curve, indicated by the light yellow area in Figure 20A. This region is determined based on vibration velocities higher than 10% of the peak vibration velocity. Notably, the extent of the pale yellow region reveals that the duration of the S-wave gradually decreases as the frequency increases, consistent with the general principle that higher frequencies decay more rapidly. However, the 30 Hz frequency exhibits an anomaly, with a significantly shorter duration compared to 40 Hz and 50 Hz. To gain further insights, Fast Fourier Transform (FFT) was conducted to analysis of the waveform within the yellow region, yielding Figure 20B. The dominant frequency of the vibration velocity at the P9 generally increases with the input frequency of the S-wave, except 30 Hz. At 30 Hz, there are four distinct peak frequencies, indicating a highly complex waveform composition compared to other frequencies. This complexity can be attributed to the fact that, as shown in Figure 15, the 30 Hz S-wave led to a substantial expansion of the roadway’s plastic zone. In essence, when the 30 Hz S-wave propagates through the roadway, it induces the most severe plastic failure, resulting in a complex waveform composition.

Figure 14 reveals an anomaly in V _pa when the roadway is exposed to P-wave activity at 60 Hz. To uncover the underlying cause of this peculiarity, we conducted a comprehensive analysis, which included detailed statistics of the plastic zones on both sides of the roadway, as well as those on the roof and floor. As depicted in Figure 21, we applied the Fish to compute the plastic zone volumes on both sides of the roadway, as well as those on the roof and floor. In Figure 21A, both V _pm and V _pn are negative at 60 Hz, and V _pm is negative at 2, 5, and 10 Hz. These negative values suggest a reduction in the plastic zone size after dynamic loading, which is quite perplexing. This phenomenon can be attributed to the stress adjustment in the zone due to dynamic loading, so that the stress conditions that would have allowed the zone to undergo plastic failure no longer exist. In addition to these anomalies, V _pm gradually increases with the increase in P-wave load frequency, while the increase in V _pn is less pronounced. Figure 21B indicates that V _sm peaks after the impact of the 30 Hz S-wave. However, when the S-wave frequency is 60 Hz, V _sm is significantly lower than at other frequencies, and V _sn is also smaller compared to other frequencies. The increase in the plastic zone induced by P-wave activity is notably less compared to that caused by S-wave activity. This discrepancy arises from the alignment of the X-axis direction of the maximum principal stress in the numerical model with the primary action direction of S-wave. In other words, S-wave enhances the maximum principal stress, while P-wave enhances the minimum principal stress. Consequently, the impact of these two waveforms on the plastic zone is distinct.

The abnormal plastic zone of the roadway induced by the 60 Hz P-wave prompted a spectrum analysis of the time-velocity curve at the P9 measurement point. This analysis is presented in Figure 22, following the same data selection principles as before. The analysis of Figure 22A reveals that the significant fluctuation time of the waveform lacks a clear pattern. The peak vibration speed at P9 tends to increase with the rising frequency, except at 70 Hz, where it decreases. In Figure 22B, the main frequency components of the waveforms at 60 Hz and 80 Hz display a distinct single peak, while those of other waveforms are more intricate, often featuring double or triple peaks. The complexity of the main frequency in the waveforms is attributed to the damage of the surrounding rock when the waveform reaches P9, leading to the generation of additional wave. Figure 14 shows that V _pa at 60 Hz and 80 Hz is notably lower compared to the surrounding frequencies, suggesting a potential correlation. It seems that when the plastic failure around the roadway becomes more severe, the detected waveform’s main frequency also becomes more intricate.

Figure 23 illustrates the maximum BSR of the roadway when subjected to S-wave and P-wave of varying frequencies. The figure demonstrates a noticeable increase in roadway rockburst potential due to dynamic loading, compared to the roadway excavation. Except for the frequencies of 10 Hz and 20Hz, the BSR induced by S-wave at other frequencies surpasses that caused by P-wave of the same frequency. This discrepancy is attributed to the primary loading direction of S-wave aligning with the maximum principal stress direction. The increase in rockburst potential associated with P-wave is relatively minor, except for frequencies of 5 Hz, 10 Hz, and 20 Hz. Interestingly, the BSR due to S-wave generally decreases with increasing frequency, except for an anomaly at 60 Hz where the BSR reaches its peak, significantly exceeding that at other adjacent frequencies. Drawing insights from Figure 21A, it becomes evident that the plastic zones induced in both the roof and floor reach their minimum extent when subjected to P-wave frequencies of 5 Hz, 10 Hz, and 20 Hz. This observation underscores a noteworthy discrepancy between the findings derived from the analysis of plastic zones and those from the BSR index. This incongruity can be attributed to that when the surrounding rock undergoes plastic deformation, its stress levels tend to decrease. Similarly, as indicated by Figure 21B, when the frequency of S-wave is 60 Hz, the increment in the plastic zones observed in the sidewalls and the top and bottom plates of the roadway is significantly less compared to the surrounding frequencies. This suggests that the 60 Hz S-wave dynamic load induces plastic deformation in a more confined region. It explains why this frequency yields the highest rockburst potential. This leads us to a plausible conclusion: the 60 Hz S-wave dynamic load elevates the maximum principal stress but does not quite reach the threshold for plastic failure in the material, resulting in a higher calculated BSR value.

Based on the engineering background of the Dongtan Coal mine, this study compared wave characteristics between near-field and far-field strong mine tremors. Additionally, investigates the propagation and attenuation characteristics of P-wave and S-wave at various frequencies, assessing their impact on the plastic failure and the rock burst risk of roadways. The main conclusions are as follows:

(1) Wave analysis reveals that compared with near-field strong mine tremors, far-field tremors are dominated by low-frequency components below 5 Hz, which may caused by the rapid attenuation of high-frequency components along with propagation.

In recent years, the frequency of far-field strong mine earthquakes has been on the rise, making the study of their impact on roadways a crucial research area. Due to the challenges in conducting accurate field or theoretical research on this issue, numerical simulation methods have become indispensable. The conclusions drawn in this study are based on the modeling of the actual engineering background conditions of Dongtan Coal Mine. It can be predicted that the simulated dynamic load propagation law and its influence on the roadway may be quite different under different conditions due to the differences in the strata structure, the properties of coal and rock mass, and the actual stress state of the roadway. Therefore, it is an effective and reliable means to establish an accurate numerical model for specific engineering research. At present, sensors for vibration velocity monitoring in engineering are arranged in the form of points, so a small amount of data collected by sensors makes it difficult to range the real situation of the site. If large-scale vibration velocity monitoring can be realized in the future, field data-based research will flourish.