1 Introduction

Located in southern Poland, the Upper Silesian Coal Basin (USCB) is an area where intensive exploitation of hard coal has been conducted for many years. This exploitation is accompanied by seismicity, induced by mining activities. This seismicity is of a twofold nature. The first type of seismic event includes typical mining seismic events, directly related to mining activities in the area of active underground workings, while the second one results from the overlapping of mining works and stresses in the tectonic active zones. This second type of mining-induced seismic event is not only particularly dangerous for people working underground but also troublesome for the inhabitants of mining areas. These high-energy seismic events have magnitudes that can exceed 4 on the Richter scale. There are relatively few high-energy seismic events among all the seismic events. In 2018, there were about 1500 seismic events with magnitude over 1.7. Among them, 29 seismic events reached the magnitude above 3. However, these are the most dangerous and cause the most damage to construction facilities.

Spectral element method (SEM) was originally developed for calculating fluid dynamics, but has also been successfully used to solve problems related to the propagation of seismic waves in rock media. SEM is mainly used in seismology for modeling wave propagation after strong earthquakes (Komatitsch et al. 2004; Maheshreddy and Raghukanth 2017). In these publications, the authors show the application of numerical methods to estimate the ground motion from earthquakes in Los Angeles Basin and Himalaya Region. Simulated ground motion is compared with the recorded data, and peak ground acceleration contour maps are provided. The simulations account for 3D variations of seismic wave velocity, density, and topography. They demonstrate that it is feasible to predict ground motion, thus improving our ability to assess seismic hazard. There are also publications presenting the results of calculations for mining seismic events in mines, e.g., zinc and copper ores (Gharti et al. 2017). In Poland, this method has not been used to study seismic effects induced by mining seismic events, but it may prove to be a very useful tool for determining the impact of vibrations caused by mining seismicity on building facilities on the surface. As a result, we obtain full waveforms, allowing effective assessment of seismic risk in the zones of mining seismic event occurrence.

The paper presents the results of modeling carried out using SpecFEM3D software (www.geodynamics.org/cig/software/specfem3d), the first authors of which were Dimitri Komatitsch and Jean-Pierre Vilotte (1995–1997), followed by Dimitri Komatitsch and Jeroen Tromp. All SpecFEM3D programs are written in Fortran 2003 (Adams et al. 2008) and use parallel programming based on the Message Passing Interface (MPI) communication protocol. The software is released under the GNU free software license (Komatitsch et al. 2018).

2 November 8, 2018, MW = 3.8 seismic event in the Budryk mine

On November 8, 2018, in the mining area of the Budryk hard coal mine located in the western part of the USCB, a seismic event with a magnitude of 3.8 on the Richter scale occurred. The geographic coordinates of the epicenter were determined on the basis of the underground seismological network and were 50.1855 N and 18.7540 E. The seismic event was located in the immediate vicinity of the Barbara fault, whose fault surface extends on the west-east axis and is inclined towards the south. It was strongly felt on the surface but the inspection of mine excavations carried out after the seismic event did not reveal any influence on the technical condition of the mining objects. After the event, about a hundred reports of possible damage to buildings were registered from residents.

For the seismic event, the tensor of the seismic moment was calculated. It describes the distribution of forces in the seismic source, as a combination of force couples with the moment of couple. Displacements in the far wave field caused by the system of forces in the source are the sum of the displacements caused by individual force couples (Aki and Richards 1980). The calculations of the seismic moment tensor were made using the program Foci (Kwiatek 2009). On the basis of the source mechanism with a relatively high shear component, it can be concluded that the cause of the seismic event was the relaxation of rock masses under tectonic stress in the Barbara fault area, to which the stresses caused by the current exploitation in this region contributed. The parameters of the mechanism of the source were determined on the basis of seismograms registered by the Upper Silesian Regional Seismological Network (www.grss.gig.eu). Registrations from three selected stations (MUR, JKW, and IGN) were used to compare the calculated synthetic seismograms with real ground motions. Figure 1 shows the area in question, along with the location of the seismic event location and the seismological sites included in the analysis.

Fig. 1
figure 1

A map of the discussed area with the location of the seismic source and surface seismic stations

Full size image

The result of the calculation of the seismic moment tensor is shown with the angular parameters of the nodal planes and compression stress (P) and tensional stress (T) in Fig. 2. This is the solution of the full tensor, for which the following are distinguished: the volume change component (ISO), the uniaxial compression or tension component (CLVD), and the shear component (DC). This seismic event was characterized by a normal slip mechanism.

Fig. 2
figure 2

Parameters for the solution of the mechanism of the focus of the seismic event (ΦA,B—azimuth of plane A,B; δA,B—plane dip A,B; λ—slip angle; ΦP,T—axis azimuth P,T; δP,T—immersion of axis P, T; ISO—percentage of the isotropic component, CLVD—percentage of the component corresponding to uniaxial compression /−/ or tension /+/, DC—percentage share of the shear component)

Full size image

The full moment tensor decomposed into the volume change component (Expl), the uniaxial compression or tension component (CLVD), and the shear component (DBCP) indicates that the solution is dominated by the shear component (62%) over the uniaxial compression (20%) and explosion (18%) components. The source mechanism points to a shear-type mechanism of detachment faulting. Nodal planes have a west-east azimuth (plane A Φ° = 269° and φ = 60° or plane B Φ° = 86 and φ = 31°). The main compressive stresses (P) operate at an angle of 75°, and tensile strains (T) operate at an angle of 14°.

3 Numerical technique of simulation

The SEM is an extension of the classical finite element method. It allows us to carry out numerical modeling of full waveforms. The SEM uses multinode finite elements, but in this method, unlike the finite element method, an unequal distribution of nodes is used in the elements, and the higher order Lagrange polynomials are used to approximate interstitial values. This approach reduces the number of nodes at the wavelength to approximately 5 and guarantees that the obtained mass matrix takes the diagonal form. The modeling of the time signals of seismic wave velocities is based on a seismic moment tensor describing the arrangement of forces in a seismic source.

The mathematical basis of the method is presented in the following papers (Komatitsch and Tromp 1999; Komatitsch and Tromp 2002a, b), and this paper presents a general outline of it. The first stage of modeling is discretization of the continuous medium to form nonoverlapping cubic elements. Each element is defined by nodes and shape functions. The shape functions approximate the displacement field within the element, whereas the system of equations that should be solved in order to determine the displacements in the nodes is calculated on the basis of the motion equation, which can be written as:

$$ Mu^{{\prime\prime} }+ Cu^{\prime }+ Ku=F $$
(1)

where u is the global displacement vector, M is the global mass matrix, C is the global edge absorption matrix (suppression), K is the global stiffness matrix, and F is the force vector.

4 Model and mesh implementation

Figure 3 shows the fragment of a grid of points constructed for calculations of waveforms as well as the values ​of density-velocity parameters for individual model layers. In the calculations, a flat-parallel velocity model with three layers was adopted. Each medium layer is described by the propagation velocities of P and S seismic waves and the density. These parameters were set on the basis of seismic measurements made by seismic tomography method in adjacent mining excavations.

Fig. 3
figure 3

Fragment (4 km × 6 km) of the grid used for numerical calculations and the density-velocity model adopted in the calculations

Full size image

The entire model covers the 49 × 49-km study area, from 49.9696 N to 50.4105 N and from 18.4571 E to 19.1425 E. The depth of the model is 2800 m. The measuring grid consists of 1,849,536 elements, each with 125 nodal points. After subtracting the common points lying on the edges of the elements, the number of nodes in the entire model is determined to be 121,857,120 (365,571,360 degrees of freedom, because the displacement components at each point of the grid are searched in three directions). The smallest and largest distances between nodes are 17 m and 51 m, respectively. Assuming that only five samples per wavelength (Seriani and Priolo 1994) are sufficient for the full waveform reproduction, we can obtain solutions for frequencies up to 5.6 Hz. In modeling, Stacey absorbing conditions defined in xmin, xmax, ymin, ymax, and zmin and free surface in zmax were applied.

The seismic source is located in the central part of the model, at a depth of 1200 m. A full moment tensor and the Gauss source function with a duration of 0.3 s were used in the calculations. The sampling rate was 625 Hz, and the length of the recording of synthetic seismograms was 35.2 s, which gives 22,000 samples for each trace.

The calculations were made on a supercomputer with four calculation nodes used for modeling. Each node was equipped with 24 processors with an allocated 5-GB memory, which in total gives 480 GB of computing memory.

5 Simulation of the November 8, 2018, MW = 3.8 seismic event

Figures 4, 5, and 6 present the results of the synthetic seismogram calculations using the SEM method in three selected points of the model, in which the correlation between recorded and synthetic seismograms was the most accurate. Next, synthetic seismograms were compared with real vibrations recorded by surface seismological stations MUR, JKW, and IGN. The epicentral distance for individual stations ranges from 20.7 to 24.7 km. They are located in different directions relative to the seismic event source, and their positions are shown in Fig. 1.

Fig. 4
figure 4

The result of modeling in the SpecFEM3D program for the MUR station (a) and the actual vibrations registered at this station (b)

Full size image
Fig. 5
figure 5

The result of modeling in the SpecFEM3D program for the JKW station (a) and the actual vibrations registered at this station (b)

Full size image
Fig. 6
figure 6

The result of modeling in the SpecFEM3D program for the IGN station (a) and the actual vibrations registered at this station (b)

Full size image

All data recorded were filtered using a Butterworth filter in the frequency range from 0 to 6 Hz. In the studied area, the terrain surface changes to a small extent, so the calculations do not take into account the influence of topography. As can be seen in the recorded data, the S wave dominates, but the P wave and surface waves are also clearly visible. The signal duration for the S wave is approximately 2 s at the JKW and IGN stations and approximately 1 s at the MUR station, which results from the smaller distance of the station from the seismic source. The peak vibration values at the JKW and IGN stations are similar and range from 1.58 to 2.28 mm/s for horizontal components and from 0.64 to 0.68 mm/s for vertical components (at epicentral distances of 22.9 km and 24.7 km, respectively). The MUR station, despite being closest (20.7 km to the east of the seismic source), is characterized by significantly lower vibration values. The horizontal components are 0.08 mm/s and 0.31 mm/s, and the vertical component is equal to 0.05 mm/s. This is due to the mechanism of the source and hence the directions of the seismic wave radiation in the rock mass.

Recordings of real vibrations have been registered on three selected measuring stations installed on the ground surface of buildings. The identification of individual wave groups is not difficult. S wave amplitudes are much higher than P wave amplitudes; in particular, at the JKW and IGN stations, the boundary between these two waves is very clear. In all records, the surface wave is also marked, although its amplitude does not exceed the amplitude of the S wave. Understandably, the calculated synthetic seismograms and the actual vibration recordings differ. This is due to many reasons. The most important reason is that the registrations took place at a relatively large distance from the source center, where ground motions due to the occurrence of wave phenomena such as dispersion or interference have a complex character. Another factor affecting the nature of the recording may be the installation of seismometers, which were installed on the buildings, and the response of the whole building to the propagating seismic wave is noted in vibration recordings. However, most importantly, the recorded vibrations are similar to the calculated synthetic seismograms, and even if the calculated signals do not fully reflect the real waveforms, they can be very useful in determining the peak vibration value ​(PGV) at any point in the model. This, in turn, can be very helpful in determining the effects of seismic-induced vibrations in areas with mining exploitation.

A quality analysis of modeling was conducted for the obtained results. It was based on the root mean square error (RMSE) criteria, which represents the square root of the average of squared differences between calculated and measured values of ground motions. The RMSE was calculated separately for P and S waves. The results of this analysis are shown in Table 1.

Table 1 RMSE of the amplitudes for P and S waves at selected stations
Full size table

Figure 7 shows the distribution of the peak resultant value of the horizontal vibration velocity component (PGVXY), and Fig. 8 shows that of the vertical component (PGVZ), calculated by the SEM method. The source is located in the center of the area shown. Visible peak values resulting from the adopted seismic focal model and in the epicenter are over 35 mm/s for the horizontal component and approximately 20 mm/s for the vertical component. As the epicentral distance increases, the peak horizontal components are located along the northwest-southeast and northeast-southwest axes and the vertical component is located along the north-south axis. In order to facilitate the verification of the obtained results, Table 2 shows recorded and calculated values of resultant horizontal and vertical components of PGV.

Fig. 7
figure 7

Distribution of the peak values of the horizontal component of the ground motion velocity PGVXY

Full size image
Fig. 8
figure 8

Distribution of the peak values of the vertical component of the ground motion velocity PGVZ

Full size image
Table 2 Comparison of recorded and modeled PGV values at selected stations
Full size table

In the SpecFEM3D program, it is not possible to take into account the amplification of vibrations by the surface layers, and this, as discussed in the literature (Olszewska and Mutke 2018), can amplify the amplitudes of vibrations several times and thus decisively influence the predicted magnitude of ground motions and the assessment of the impact of seismicity on buildings at ground level. The value of the vibration amplification coefficient can be determined in several ways: by using an analytical approach, by direct measurement of the amplification effect, or by testing the spectral properties of the recorded seismometer signals. The most commonly used method in the USCB is the analytical method proposed by Savarienskij (1959). The model adopted in this method is a 1D model, and the knowledge of such parameters of the medium, as the longitudinal wave velocity in the rock bed and overburden, the density of the substrate and overburden, the thickness of the overburden, and the damping factor, is necessary for the calculation. As a result, we obtain the distribution of the vibration amplification coefficient as a function of frequency. By knowing the dominant frequencies of vibrations reaching the overburden layers, maps of the distribution of the amplification coefficient can be created. Then, by multiplying the amplification coefficient and vibration values calculated by the numerical methods, the predicted vibrations can finally be obtained throughout the studied area.

To describe the effects of mining seismic events on the surface, the Mining Seismic Instrumental Intensity Scale (MSIIS-15) is applied in the USCB. The MSIIS combines instrumentally measured parameters of vibrations with observed macroseismic effects. These macroseismic observations are arranged in six descriptive seismic intensity degrees. MSIIS-15 is a two-parameter scale, which assigns seismic events induced by mining to certain intensity degrees. The assignment is based on the peak amplitude of horizontal vibration velocity and the duration time of signal. The MSIIS-15 allows an approximate assessment of the impact of mining seismic events on buildings, within the range from harmless vibrations up to slight structural damages (Mutke et al. 2015). The distribution of individual degrees of the MSIIS-15 is shown in Fig. 9.

Fig. 9
figure 9

The Mining Seismic Instrumental Intensity Scale (MSIIS-15)

Full size image

For all three stations analyzed in the article, the PGV values of horizontal components are below 5 mm/s. That means that the effects of vibration are within the first degree of the MSIIS-15. This is due to the long epicentral distance. But it is possible to show the distribution of intensity degrees for the whole area of interest. Such a map for the Budryk mine area is shown in Fig. 10. The map was made assuming signal durations above 3 s (then, we are dealing with a higher harmfulness of vibrations).

Fig. 10
figure 10

Distribution map of intensity for November 8, 2018, seismic event

Full size image

The MSIIS-15 is very useful. It allows us to assess the predicted damage distribution on the surface according to numerous seismic observations in the mining regions. This is why this scale is of great practical significance, because it helps estimate the seismic hazard in mining areas.

6 Conclusions

This paper presents the possibilities of using numerical modeling with the SEM to determine ground motion parameters from seismic events induced by mining activities. For the USCB area, such modeling was carried out for the first time. The analytical methods used thus far were based on the determination of the relationship (damping relationship) between the parameters of ground motion velocity or acceleration, and factors affecting the size of these parameters, such as seismic event energy, distance from the source, or soil class. The use of numerical modeling allows us to include a complex mechanism of the seismic source in the calculations by applying the tensor of the seismic moment, which best reflects the distribution of forces in the source. Synthetic seismograms are calculated based on wave equations and enable full waveform imaging, while damping relationships are based on a statistical approach. The obtained synthetic seismograms do not fully reflect the wave motion recorded by surface seismological stations. This particularly applies to records from stations located at a considerable distance from the seismic source, where ground motions have a complex character. The peak ground velocity maps presented in the article do not include the topography and they were generated assuming the homogenous ground surface. Topography is not an obstacle in this case, because a flat topography is concerned, but the assumption of homogenous ground surface is a significant simplification. The thickness and physical parameters of the surface layers for the whole model would have to be known to take into account a near-surface geological building. However, synthetic seismograms provide the opportunity to determine the peak vibration values and, as a result, they can complement analytical methods in solving issues related to the prediction of seismic effects in areas with mining activity. Synthetic seismograms can be helpful in improving the efficiency of such a forecast and thus contribute to improving assessment of seismic hazard by referencing previous seismic events, as well as the effects of planned exploitation.