Induced Seismic Events—Distribution of Ground Surface Displacements Based on InSAR Methods and Mogi and Yang Models

Abstract

In this article, we present a possible approach to use satellite radar data for a complete description of the formation process of a subsidence trough resulting from an induced seismic event—a mining tremor. Our main goal was to verify whether SAR data allow for the calculation of the basic indicators for the trough (w—subsidence, T—trough slope, K—curvature, u—horizontal displacements, $$\epsilon$$—horizontal deformations). We verified the extent to which the Mogi and Yang models can be fitted to match the actual displacements recorded after an induced seismic tremor. The calculations were performed for the Legnica-Glogow Copper Belt (LGCB) area in southwest Poland. Due to intensive mining operations and specific geological and tectonic conditions, the area shows a high level of induced seismic activity. Our detailed analysis focused on four powerful mining tremors: the first tremor occurred on 29 November 2016 (M$${}_W$$3.4), the second on 7 December 2017 (M$${}_W$$3.3), the next on 26 December 2017 (M$${}_W$$3.6) and the last tremor on 29 January 2019 (M$${}_W$$3.7). For each analyzed event, we determined the displacements based on the Differential Interferometric Synthetic Aperture Radar (DInSAR) method and Sentinel 1 synthetic aperture radar (SAR) data from two paths (22 and 73). Additionally, for the period from November 2014 to October 2020, we calculated the displacements using the Small Baseline Subset method (SBAS) time series method. In all cases, the tremor was followed by the development of long-lasting surface deformations. The obtained results allowed us to conclude that it is possible to calculate indicators that result from a specific induced mining event. Considering the full moment tensor and nature of the tremor source, we demonstrated that the Mogi and Yang models can be employed to describe the influence of an induced tremor on the surface in an area of mining activity. We also confirmed the global character of the influence of the reduced troposphere on SAR data calculations. Our conclusions indicate that accounting for the tropospheric correction does not distort horizontal and vertical displacement values in regions influenced by mining activity/tremors.

1. Introduction

Examples provided to date in the literature indicate that synthetic aperture radar interferometry (InSAR) allows for the detection of local ground surface changes. Based on the SAR satellites that are currently active, we can calculate displacements of the Earth’s crust caused by such phenomena as mining activity (both underground and on the surface) [1,2], landslides [3,4,5], oil and gas extraction [6,7], pumping water from underground reservoirs [8,9] and floods [10,11]. These phenomena also include induced seismic events, which—if powerful enough—may lead to the development of long-lasting displacements on the ground surface. In all these cases, the induced seismic events are directly caused by human manipulation of the rock mass (see Table 1, the “type of event” column). These events result from an imbalance in the rock mass, which leads to the release of accumulated potential energy. Some small part of this energy is transformed into seismic energy, which propagates from the tremor center in the form of elastic waves. In contrast to natural seismic events, source of induced tremors have an isotropic (non-double couple) character and a much lower force and affect areas which are several times smaller.

To date, research on the application of SAR data to induced seismicity has focused mainly on the calculation of Line-of-Sight (LOS) displacements [15,26,27], and more sporadically on the calculations of vertical and horizontal displacements [28,29]. Our main goal was to verify whether synthetic aperture radar interferometry data allow for the calculations of the basic indicators (w—subsidence, T—trough slope, K—curvature, u—horizontal displacements, e—horizontal deformations) that describe a subsidence trough that has formed due to a mining tremor. We also wanted to demonstrate that the use of tropospheric correction has an influence on the quality of the calculated displacement components.

In this article, we focused on analyzing induced mining tremors that were recorded in 2016, 2017 and 2019 in the area of the Legnica-Glogow Copper Belt (LGCB) in Poland. Our primary goal was to use radar data from Sentinel 1A/B satellites to calculate LOS displacements and to then calculate vertical displacements, as well as displacements in the E-W direction. We compared our results with results of the calculations based on the Mogi and Yang theoretical models. We further calculated the basic indicators ($w,T,K,u,e$) for selected cross-sections of the formed subsidence troughs, as these serve to define the categories of mining areas.

Milczarek [30] demonstrated that it is possible to use the Mogi model for the modeling of vertical displacements due to mining activity on the ground surface. In this article, we expand our research to the use of Mogi and Yang models and satellite radar interferometry to determine the displacements (vertical and horizontal) caused by a mining tremor. By definition, both models are used to describe the magnitude and the direction of ground displacements due to changed volume or the pressure of a point in a uniform elastic half space. Referring to our own calculations, we believe that the models can be employed to model the influence of induced mining tremors on the ground surface.

2. Background and Methods

Ground surface displacements caused by a seismic event (natural or induced) can be calculated, e.g., on the basis of analytical models. In areas with no tectonic activity, mining-induced seismic events cause ground surface displacements, the distribution of which fits the influence function theory by showing a normal distribution [31,32]. In our previous articles [30], we demonstrated that deformations related to mining-induced tremors are located in the already existing terrain subsidence areas. These subsidences are in turn the result of underground mining activity. Therefore, tremors observed in the area of mining operations intensify the anthropogenic impact on the ground surface. We believe that, in some cases, simulations of tremor-related displacements in mining areas can be carried out with the use of models that are applied in the case of volcanic activity.

2.1. Mining-Induced Seismicity—Source Mechanism

The nature of the earthquake source mechanism is determined as a moment tensor, which can be represented by matrix M:

(1)$M=\left[\begin{array}{ccc}{M}_{11}& M_{12}& M_{13}\\ M_{21}& M_{22}& M_{23}\\ M_{31}& M_{32}& M_{33}\end{array}\right]$

The moment tensor is decomposed into three components: $M_1$, which determines the volume changes in the source; $M_2$, which is responsible for the source rotation of the stationary medium; and $M_3$, which can be divided into a compensated linear vector dipole (CLVD) and a shifting movement caused by the double couple forces (DC). The seismic moment tensor is determined by inversion based on recorded seismograms. The inversion results are presented in three forms: a full moment tensor, including volume change (implosion, explosion); the CLVD component, which is responsible for uniaxial compression or extension; and the shear component, which is described by the double couple forces. The percentage share of the individual components is also determined. This fact allows us to identify which process dominates in the source [33]. Sen et al. [34] state that most natural earthquakes do not demonstrate isotropic or uniaxial compression. However, in the case of mining-induced seismic events, the non-double couple components are often significant or dominant. Non-double couple components are related to such phenomena as tremors resulting from the caving in of a ceiling over the extraction area or from the deflection and cracking of the ceiling and tremors associated with the cracking of the protective pillar. Determining the mechanism of the tremor source induced by mining operations is significant in the case of modeling the impact of these phenomena on the terrain surface.

2.2. The Mogi Model

The first model used to verify the results of InSAR calculations was the Mogi model [35], which was originally developed to describe deformations caused by volcanic activity. However, the universal approach allows the estimation of ground displacements resulting from various phenomena that are related not only to volcanic activity [36,37] but also to underground gas extraction and storage [38] and to metal ore mining [30]. The Mogi model is based on two basic assumptions [35]: the Earth’s crust is a partially infinite elastic body, and its surface is subject to displacement due to the changing pressure or volume of the spherical source (e.g., a magma chamber or a mining excavation) (Figure 1).

2.3. The Yang Model (Prolate Spheroid Model)

The second of the models used in the verification of the obtained calculation results was the Yang model [39], which is a development on the idea of a symmetrical, spherical point source implemented in the Mogi model. It was introduced with the aim of better fitting the source model to the borders of the magma chamber. Therefore, it is most typically used in the calculations of volcanic activity [40,41] but can also be applied in such tasks as calculating displacements in the production fields of geothermal power plants [42,43]. The following assumptions were made in the Yang model [39]: the displacement of an elastic half space—i.e., the Earth’s crust—is due to the change of pressure exerted on an elongated spheroidal source that is submerged in it and which has finite borders (Figure 1). The identified displacement is the sum of two components, of which the first is related to the position of the center of body dilatation and the second is related to double forces that generate uniform normal stresses on the surface of the spheroid.

2.4. DInSAR—Horizontal Displacements

Ground surface displacements caused by natural and induced phenomena are characteristic in that—depending on the type of the phenomenon—the displacement may be observed in both the vertical and the horizontal plane. For this reason, in the case of areas affected by earthquakes or by underground mining activity, it is important to identify the extent of displacement in both planes.

DInSAR-calculated ground surface displacements are in LOS. The LOS displacement vector deviates from the vertical by the incidence angle $\theta$. The value of the incidence angle varies for individual SAR sensors. For an individual Sentinel-1 image, the value may be from approximately 29${}^{\circ }$ for near-range pixels up to 46${}^{\circ }$ for far-range pixels. Furthermore, the trajectory of the orbiting satellite is off the northern direction by an azimuth $\alpha$ (heading angle). In the case of the Sentinel-1 satellite, this value is approximately −15${}^{\circ }$ for the ascending orbit, and approximately −165${}^{\circ }$ for the descending orbit. The geometry of the SAR-based LOS displacement measurement from two orbits is shown in Figure 2.

The LOS displacement vector $d_{LOS}$ comprises three components: the vertical $d_V$, the north–south horizontal (N-S) $d_{NS}$ and the east–west horizontal (E-W) $d_{EW}$. The relationship between the vector and the three components can be represented as follows:

(2)$d_{LOS}=\left[\begin{array}{ccc}-sin\theta cos\alpha & sin\theta sin\alpha & cos\theta \end{array}\right]\left[\begin{array}{c}{d}_{EW}\\ d_{NS}\\ d_V\end{array}\right]$

With the displacement data available from several sources—e.g., from various SAR sensors or from several independent geometries—the above equation can be represented as follows:

(3)$\stackrel{\mathbf{v}}{\left[\begin{array}{c}{d}_{LOS1}\\ d_{LOS2}\\ ⋮\\ d_{LOSn}\end{array}\right]}=\stackrel{\mathbf{A}}{\left[\begin{array}{ccc}-sin{\theta }_{1}cos{\alpha }_1& sin{\theta }_{1}sin{\alpha }_1& cos{\theta }_1\\ -sin{\theta }_{2}cos{\alpha }_2& sin{\theta }_{2}sin{\alpha }_2& cos{\theta }_2\\ ⋮& ⋮& ⋮\\ -sin{\theta }_{n}cos{\alpha }_n& sin{\theta }_{n}sin{\alpha }_n& cos{\theta }_n\end{array}\right]}\stackrel{\mathbf{x}}{\left[\begin{array}{c}{d}_{EW}\\ d_{NS}\\ d_V\end{array}\right]}$

The above equation indicates that a minimum of three independent data sources are needed in order to calculate the displacement components in all three directions. Because of the heading angle $\alpha$ values assumed for Sentinel-1, the N-S component is estimated with a significant error, even if data from several different geometries are available. In the case of Sentinel-1, a common approach is to ignore the N-S component when deriving the displacement components from multiple geometries [44].

After ignoring the N-S component, the vertical and the E-W components can be calculated with the use of LOS displacements from two paths: one ascending and one descending. The second column of matrix A in Equation (3), which describes the N-S displacement, can thus be removed. In this case, the displacements from two paths may be represented for each pixel with the following equation:

(4)$\left[\begin{array}{c}{d}_{asc}\\ d_{dsc}\end{array}\right]=\left[\begin{array}{cc}-sin{\theta }_{asc}cos{\alpha }_{asc}& cos{\theta }_{asc}\\ -sin{\theta }_{dsc}cos{\alpha }_{dsc}& cos{\theta }_{dsc}\end{array}\right]\left[\begin{array}{c}{d}_{EW}\\ d_V\end{array}\right]$

For each pixel, the $\theta$ incidence angle in the equation above should assume certain values resulting from the imaging geometry; i.e., lower values for the near-range pixels and higher values for the far-range pixels. Displacement components are obtained by inverting the above equation.

For most areas around the globe, the Sentinel-1 system provides coverage with at least two images: one from the ascending path and one from the descending path. With the use of these data, we can thus calculate the vertical and the E-W horizontal displacement components for a particular area.

In the case of displacements resulting from a natural earthquake, which can be typically classified as non-isotropic, it is necessary to identify all three displacement components. An approximate estimation of the N-S component is possible among others with the use of measurement results obtained with the azimuth pixel offset method. If the Pixel Offset results from the ascending path $d_{PO-asc}$ and from the descending path $d_{PO-dsc}$ are included in Equation (4), it will be extended as follows [45]:

(5)$\left[\begin{array}{c}{d}_{asc}\\ d_{asc}\\ d_{PO-asc}\\ d_{PO-dsc}\end{array}\right]=\left[\begin{array}{ccc}-sin{\theta }_{asc}cos{\alpha }_{asc}& sin{\theta }_{asc}sin{\alpha }_{asc}& cos{\theta }_{asc}\\ -sin{\theta }_{dsc}cos{\alpha }_{dsc}& sin{\theta }_{dsc}sin{\alpha }_{dsc}& cos{\theta }_{dsc}\\ sin{\alpha }_{asc}& cos{\alpha }_{asc}& 0\\ sin{\alpha }_{dsc}& cos{\alpha }_{dsc}& 0\end{array}\right]\left[\begin{array}{c}{d}_{EW}\\ d_{NS}\\ d_V\end{array}\right]$

2.5. Tropospheric Delay Effect

The troposphere, the lowest and the most dynamic layer of the Earth’s atmosphere, has a significant influence on the passing SAR signal. Tropospheric refraction is determined by the spatial distribution of temperature, pressure and air humidity. It can be classified into two components: the dry (hydrostatic) component and the wet component [46]:

(6)$N={\left(k_1\frac{P}{T}\right)}_{hydr}+{\left(k_2\frac{e}{T}+k_3\frac{e}{T^2}\right)}_{wet}=N_{hydr}+N_{wet}$

where P is dry air pressure (hPa); T is temperature (K); e is water vapor pressure (hPa); and coefficients $k_1$, $k_2$, $k_3$ are empirically determined constants, whose values are, respectively, $k_1$ = 77.6 K · hPa${}^{-1}$, $k_2$ = 23.3 K · hPa${}^{-1}$ and $k_3$ = 3.75 K${}^2$· hPa${}^{-1}$ [47].

Due to tropospheric refraction, the signal phase is delayed. The phase component of the signal representing the tropospheric delay (${\varphi }_{tropo}$) is described by the formula:

(7)${\varphi }_{tropo}=\frac{4\pi }{\lambda }·\Delta R_{tropo}$

where $\lambda$ is the electromagnetic wavelength and $\Delta R_{tropo}=R_1-R_2$ is the change in the range of the electromagnetic wave due to the tropospheric delay between two image acquisitions.

Two types of delay are identified: turbulent delay due to such phenomena as thermal convection and differences in wind speed and wind direction on different altitudes, friction and complex weather patterns; and stratigraphic delay, correlated with the terrain topography [48]. Various approaches allow the influence of tropospheric disturbances on InSAR measurements to be reduced: empirical methods [49,50], methods based on data from weather models [51,52], spectrometry methods [53] and methods employing Global Navigation Satellite System (GNSS) measurements [54].

In this case, the influence of the troposphere was reduced with the use of the Generic Atmospheric Correction Model (GACOS) [55], which is based on the Iterative Tropospheric Decomposition Model (ITD) [56]. In this approach, the delay is calculated iteratively by decomposing the signal into stratigraphic and turbulent components.

The tropospheric correction is allowed for when calculating LOS displacements, which in turn serve to calculate the vertical and the E-W horizontal displacement components.

2.6. Subsidence Trough Indicators

Ground surface deformations caused by mining activity manifest on the surface in the form of subsidence troughs. The parameters of the trough are defined and characterized by a number of geometric values; i.e., the deformation indicators described in Table 2. Figure 3 graphically shows the relationships between individual indicators. These deformation indicators serve, among other purposes, to determine the categories of mining areas, which in turn describe the influence of mining activity on the ground surface and on the infrastructure. The indicators may also serve to provide a more detailed description of changes that occur within a subsidence trough over time and/or of the influence of induced seismicity-related phenomena on the development of the trough.

Ordinarily, the indices are determined on the basis of cyclic, classical geodetic (in situ) measurements of leveling and tachymetric surveying, which are used to determine changes in the height of measurement points and changes in the length of sides along the observation lines (Table 3). In our approach, we show that the determination of indices, concerning short-term subsidence and horizontal displacements (in the E-W plane), is possible using DInSAR measurement data and LOS vector decomposition. The subsidence and horizontal displacements are determined directly from the obtained measurement results, while the slope and strain are determined as the change in subsidence and horizontal displacement for a particular distance.

3. Application Examples

Seismic events in the study area are induced by mining activity. The scheme of calculating SAR data is presented in (Figure 4). Our analysis of induced seismic events was carried out on the basis of DInSAR and SBAS calculations for a region that undergoes underground copper ore extraction in southwest Poland (Figure 5).

We analyzed four strong mining tremors that occurred in 2016, 2017 and 2019. In each of the cases, the events had enough force to cause long-lasting displacements on the ground surface. In our previous publications, we indicated that, in this region, long-lasting ground surface displacements are observed in the case of events of at least M${}_W$3.7 [30].

In the final calculation stage, after phase unwrapping, we allowed for the influence of the troposphere by using the GACOS data [55]. Table 3 comprises basic information regarding the used SAR images.

Legnica-Glogow Copper Belt

The Legnica-Glogow Copper Belt is located in the area of the Fore-Sudetic Monocline in southwest Poland. The area is rich in copper ore, which is currently mined from six deposits by three mining plants—Lubin, Rudna and Polkowice-Sieroszowice—at depths which occasionally exceed 1200 m. The mining operations are performed with the retreat room and pillar method with hydraulic backfill (in the case of areas requiring surface protection). Copper ore in the LGCB region is generally formed in three types of rocks: Rotliegend and Grauliegend sandstones, Zechstein copper shales and carbonate rocks. Rocks of the last type, which are typically found in the roofs of the excavations, show high strength and an ability to accumulate elastic energy, and they therefore facilitate the occurrence of tremors.

Induced seismicity in the LGCB region is influenced by a number of factors, which include

• Complicated geological and mining conditions, such as fault and fold structures or tectonically disturbed surfaces;

• The deposit depth and its thickness;

• The presence of rocks prone to violent energy releases;

• The intensiveness of drilling and blasting operations using specialized heavy duty machinery.

Induced seismicity itself is not responsible for the formation of new subsidence troughs. However, it frequently causes significant subsidence in already existent troughs.

Our detailed analysis focused on four powerful mining tremors: the first tremor on 29 November 2016 (M${}_W$3.4), the second on 7 December 2017 (M${}_W$3.3), the next on 26 December 2017 (M${}_W$3.6), and the last tremor on 29 January 2019 (M${}_W$3.7). In each case, the tremor was followed by the development of long-lasting surface deformations. Figure 7 shows a comparison of the DInSAR-based results with allowance made for the influence of the troposphere.

The center of the trough which subsided due to the tremor on 29 November 2016 is located about one kilometer west of the western side of the Zelazny Most tailings pond. In the image of the vertical component, the trough assumes the shape of an ellipse elongated in the north–south direction. The calculated component indicates the maximum subsidence value in the analyzed period at −82 mm in the vicinity of the center of the trough. The value is greater than the value read from the LOS results for the ascending path (−73 mm) and for the descending path (−67 mm). The image of the E-W horizontal component indicates horizontal displacements within the subsidence trough at +/−30 mm over the investigated period. The ground surface was displaced towards the center of the subsidence trough. In the eastern part of the trough, the displacement in the western direction reached the western wall of the tailings pond.

The trough affected by the displacements resulting from the tremor on 7 December 2017 is located about 4 km west of the trough analyzed in the case of the 2016 tremor. The effects of this tremor, as in other cases, were observed within the existing subsidence trough. The maximum subsidence value in the analyzed period was at −75 mm, with the LOS values being −60 mm ascending and −67 mm descending.

The image of the vertical displacement component indicates that the trough is circular in shape, with the maximum subsidence value in the center recorded at −83 mm (ascending: −70 mm, descending: −88 mm). The E-W horizontal component suggests, as in the case of the 2016 tremor, that the ground became displaced within the trough towards its center by +/−40 mm. In the case of this tremor, the non-uniform horizontal displacement distribution and the varying subsidence values obtained from the two paths may suggest that the displacement within the trough was less uniform than in the case of the 2016 tremor.

The effects of the tremor on 29 January 2019 were observed approximately 2 km east of the subsidence troughs caused by one of the 2017 tremors. As in the two previous cases, the mining tremor occurred in the Rudna Glowna fault (NW-SE) and the Biedrzychowa fault (SW-NE) system. This region is among the most difficult in the LGCB and is characterized by very high anthropogenic seismic activity. In the image of the vertical component, the trough assumes the shape of an ellipse elongated in the N-S direction, with a maximum subsidence at the center equal to −112 mm. The horizontal components in the E-W direction range from −63 mm (east side) to +75 mm (west side).

Our current and previous research [26,30,30] shows that all troughs in the LGCB region are characterized by symmetry in at least one direction (shape of a circle or ellipse). Surface displacements caused by induced seismic events occur only within pre-existing subsidence troughs.

The investigation of the nature of seismic phenomena in the LGCB area indicates in many cases, especially those with powerful events, that the source has an isotropic character [33,58,59]. We calculated the full seismic moment tensor for all analyzed events to determine the nature of their sources. We used data registered by the local LUMINEOS seismological network [60]. All analyzed shocks are partly isotropic; thus, we decided that modeling terrain surface displacement using models describing the isotropic nature of the source was justified in these cases. Table 4 describes the characteristics of tremors based on the moment tensor, the percentage values of individual components (ISO, CLVD and DC) and figures displaying focal mechanism.

The explosive or implosive model of the study corresponds to the process of the volumetric destruction of the structure of the rock mass. This may be a consequence of the implementation of active rockburst prevention (blasting) or may possibly occur as a result of pressure on the overlying rock layers on the bed. The mechanism described by the double couple forces (DC) corresponds to the shocks associated with the cracking of thick, compact rock complexes that have high stiffness and strength and which occur in the roof of the bed, as well as the movement of rock masses on the fault planes. The analysis of the components of the seismic tensor (Table 4) showed that, in the case of the analyzed tremors, the explosive component ranged from 26 to 33%, and the CLVD ranged from 52 to 65%. The share of the double couple forces was significantly lower, from 9 to 22%. Such a distribution of tensor components applies to volumetric (V) events and to pillar cracking (CLVD) events, which are simultaneous.

The next step was to find a fit between an appropriate model and the results obtained with the DInSAR method. The initial fitting of the model to the selected event was based on the analysis of the shape of the calculated vertical displacements. The Mogi model was used when the dimensions of the displacement range were close to a circle. In contrast, the Yang model was chosen when the shape of the influences was elliptical. Table 5 shows the source parameters for each analyzed event, determined from DInSAR results.

The area affected on the surface by the tremor of 26 December 2017 was similar to a circle (Figure 8). Therefore, the event best fitted the Mogi model. In other cases, a better fit was obtained with the use of the Yang model (Figure 8). This was because, instead of being symmetrical, the influence of this tremor on the surface was elliptic (with a certain rotation angle). For this reason, the Yang model, which is based on an elongated spheroidal source, seemed to be the optimal choice.

The final result of our research was to find the basic momentary/temporary indicators describing the displacement caused by an induced seismic event in the region of the subsidence trough (see Section 2.6). Based on the described indicators, a category of a mining area is defined. This category is in turn a measure of the hazard due to the observed influence of mining activity on the ground surface and on the buildings. We compared the indicators calculated on the basis of the radar data with the indicators calculated on the basis of the fit models: Mogi and Yang.

4. Discussion

Due to the fact that the shocks are directly related to exploited copper ore, it was decided to carry out additional calculations using the SBAS time series method. The calculation period covered a period from 15 November 2014 to 26 October 2020. A total of 296 SAR images from track no. 73 were used for the calculations. The average period between each acquisition was 6 days, excluding the time when only Sentinel 1A orbited. The total number of interferograms was 1147. A preliminary analysis of the impact of the tremors on the surface demonstrated that the ground subsidence in the LGCB region occurred within the borders of the already existing subsidence troughs (Figure 9). The location of the epicenter also indicated a correlation with the performed mining works (Figure 5).

In the LGCB region, copper ore exploitation from a single bed and with the use of the room and pillar method causes an approximately time-constant ground surface subsidence. In such cases, the influence of induced seismic events can be easily observed (Figure 9). The increments of ground subsidence resulting from such events do not affect the shape of the trough. Instead, their isolated occurrences result in rapid trough developments. By allowing for the tropospheric component in the calculations, it was possible to significantly improve the results. The clear influence of the tropospheric delay has a global character (Figure 10). Allowance for this delay reduces the displacements within the entire radar image by a certain constant value.

One of the key elements of our research was to match an appropriate analytical model to the observed influence of selected mining tremors on the ground surface (Figure 10). In Section 3, we present the results of fitting the Yang model (for the tremor events from 7 December 2016, 7 December 2017 and 29 January 2019) and the Mogi model (for the 26 December 2017 event). We believe that the calculated residua for the vertical and the horizontal displacements demonstrated that the models were fitted properly to each of the analyzed events. In all cases, we assumed the location of the sources (spheroidal and elliptic) to be at a depth corresponding to the actual mining depth. Nevertheless, it should be noted that for the horizontal displacements observed in the basins, the horizontal displacement vector should be directed towards the center of the basin. For this reason, the horizontal component (EW) only represents the actual displacement in the place of the EW axis, which passes simultaneously through the center of the trough. At all points north and south of the EW axis, the horizontal displacement values are the component of the total horizontal displacement. In the case of all the tremors, the interferograms clearly show an area east of the troughs (Figure 9), which is the Zelazny Most tailings pond—the largest engineered structure of this type in Europe. The dynamics of the area are related to its function. Constant changes of the tailings, as well as works performed on the embankments, are reflected in the image of the phase of the interferometry wave. Because of the limitations of InSAR methods, the displacements observed in the region are definitely not reliable.

The displacements after the 2016 tremor were located very close (about 300 m) to the western wall of the tailings pond. In our opinion, the Yang model seems to offer a better fit in the mining activity region. This is clear in the case of the horizontal displacements (see the residua of the horizontal displacements in Figure 10 and Figure 11). Theoretically, the Yang model is also more universal; it can be successfully used for local troughs that have an elongated, elliptic shape, as well as for troughs that have a near-circular shape. However, it should be noted that the models we decided to use have their limitations, which can be significant in certain conditions. Firstly, in the case of dipping beds, it is impossible to allow for the source, which would be inclined in the vertical plane (this applies to the Yang model). The influence of mining activity in a dipping bed region is different than the influence in a region with flat strata. In such cases, the curvature of the subsidence trough is shifted by a certain value: the so-called deviation angle. The second major limitation is that tectonic conditions in the region of mining activity cannot be accounted for. Local tectonic structures may limit or amplify the influence of mining activity on the ground surface. This is of special importance with regards to the analyzed area, as the seismicity of the LGCB region proves to be directly linked to the tectonics. The main tectonic elements in the region are the Biedrzychowa fault and the Rudna Glowna fault (they have a character of a dislocation with a significant share of the horizontal component) (Figure 9).

Theoretically, the values of ground surface indicators depend on the distance from the center of the subsidence trough. This is valid, except for vertical displacements, which are the largest over the central part of the area of mining operations. The rest of the indicators reach extreme values in the vicinity of the border. The steepest slopes are at the same time observed directly over the border, and the curvatures, displacements and horizontal deformations are located at a distance from it.

Figure 11 shows the values of the indicators calculated from the SAR data, as well as the indicators calculated from the fitted analytical models. We observed a high correlation between the curves representing the vertical displacements for all the analyzed events. In the case of the 26 December 2017 tremor, the curve and the $w_{max}$ from the SAR data and from the model are very similar. In the case of the 2016 tremor, the actual $w_{max}$ is shifted to the east by approximately 200 m. Apart from the shifted $w_{max}$, the vertical displacement curves are practically identical. A better fit for the model was impossible, as the trough has an asymmetrical cross-section in the E-W direction.

The slopes and curvatures of the trough represent changing subsidence, while the deformations reflect changing horizontal displacements. The plots of the slopes, curvatures and deformations for the InSAR data show great variability. However, the main trend in the changing values is consistent with the results obtained from the theoretical model and with the principles of the influence function theory. Indicators plotted for the data from the model always have regular curves, as the model does not reflect ground surface changes with sufficient accuracy. Plots based on classic geodesic measurements do not show such variability. This is because survey data are collected with a lower frequency. Rapidly changing and/or high local values of slopes, curvatures and deformations may represent very important information. Such information is key in the identification of areas in which buildings or infrastructure may be at a particular risk of damage. At the same time, the possibility should not be dismissed that some subsidence and horizontal displacement values calculated on the basis of the InSAR methods include errors due to various decorrelations. In this case, the indicators that describe the changes of a particular parameter (slope, curvature and deformation) may prove helpful in identifying such areas in the radar interferogram. In order to identify the source of random dispersion, it is necessary to observe the subsidence trough and the distribution of indicators on a long-term basis.

5. Conclusions

In this article, we present a possible approach to using satellite radar data for a complete description of the formation process of a subsidence trough resulting from an induced seismic event—a mining tremor. We demonstrated the possibility of using SAR data to calculate the basic indicators that describe a subsidence trough. Our second objective was to verify whether the Mogi and Yang models can be employed to describe the influence of an induced tremor on the ground surface in an area of mining activity. Our third objective was to verify the extent to which allowing for the tropospheric correction in the calculations influences the values of the discussed indicators. Our research area was a region of underground copper ore extraction (the LGCB) located in southwest Poland. The region is characteristic in that it shows a high level of induced seismic activity resulting directly from underground mining operations and from geological conditions.

Based on the obtained results and on the above discussion, the authors propose several basic conclusions.

• Due to the relatively small sizes of the subsidence troughs, it was difficult to calculate all the displacement components (N-S, E-W, vertical). In most regions around the globe, SAR data coverage from Sentinel 1A/B satellites is from one or two paths (ascending and descending), and on far fewer occasions from three paths or more. The lack of data from other sources significantly obstructs the calculation of all of the components. As the tremors were local (the trough range was several hundred meters) the Azimuth Pixel Offset or Multiple Aperture Interferometry (MAI) methods did not provide satisfying results.

• The global character of reducing the influence of the troposphere significantly improves the displacement results. Moreover, it does not distort horizontal and vertical displacement values in regions influenced by mining activity/tremors.

• The presented analysis of the tremors demonstrates that they do not generate new subsidence troughs and only rapidly influence the development of already existing troughs.

• We demonstrated that if the source of an induced tremor has an isotropic character, it is possible to effectively model terrain displacements using simple isotropic models such as Mogi and Yang.

• We believe that it is relatively easy to fit the analytical models (Mogi and Yang) to the displacement image resulting from a mining tremor. We can use this as a basis for calculating the theoretical indicators describing the subsidence trough.

• A comparison of the indicators from the models with the actual data allows for a more accurate fitting of the theoretical model. The deformation indicators calculated from DInSAR measurements allow for a detailed observation of the changes that occur on the border of the subsidence trough and that are caused by mining activity and by phenomena related to induced seismicity.

• Determination of short-term trough deformation indices from DInSAR measurements and LOS decomposition is possible. Additionally, due to their large measurement range, indicators can be determined for many profiles that are located in different parts of the subsidence basin. Most importantly, the measurement is remote and does not require expensive and lengthy in-situ works.

Author Contributions

Conceptualization, W.M., A.K. and D.G.; Formal analysis, W.M., A.K., D.G. and N.B.; Investigation, W.M., A.K., D.G. and N.B.; Methodology, W.M., A.K. and D.G.; Project administration, W.M.; Resources, W.M.; Supervision, W.M.; Validation, W.M., A.K. and D.G.; Visualization, W.M., A.K., D.G. and N.B.; Writing—original draft, W.M., A.K., D.G. and N.B.; Writing—review and editing, W.M., A.K., D.G. and N.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The Sentinel-1 data presented in this study are openly available in the Earth Observing System Data and Information System (EOSDIS) (https://search.asf.alaska.edu/ accessed on 3 November 2020). In turn, the seismic data are from the open access databases: the European Plate Observing System (EPOS) (https://tcs.ah-epos.eu/ accessed on 9 November 2020).

Acknowledgments

Calculations have been carried out using resources provided by Wroclaw Centre for Networking and Supercomputing (http://wcss.pl), grant No. 345.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Schematic diagrams of the Mogi (left) and the Yang (right) models.

Figure 2. Geometry of a theoretical displacement vector dTRUE as seen by a SAR satellite in an ascending and descending pass, respectively.

Figure 3. Deformation indicators describing ground movement over an area of underground mining operations (grey area). Left-hand side: vertical components; right-hand side: horizontal components.

Figure 4. Calculation schema of the basic indicators for the subsidence trough based on SAR data with tropospheric correction (Generic Atmospheric Correction Model (GACOS)) and theoretical Mogi and Yang models.

Figure 5. Ground coverage of the Sentinel 1A/B acquisitions for the Legnica-Glogow Copper Belt area and the location of the seismic network stations (left). On the right: location of over 6000 mining tremors (MW between 0.9 and 4.0) recorded in the period from 10 December 2013 to 31 August 2020 (source: IS EPOS (2017), episode: LGCD, https://tcs.ah-epos.eu/#episode:LGCD (accessed on 9 November 2020), doi:10.25171/InstGeoph_PAS_ISEPOS-2017-006). The dashed green line represents the cross-section of the profile shown in Figure 6.

Figure 6. Location (depth) of the induced seismic events along the selected profile (Figure 5): the red line represents the exploited copper deposit (left). On the right: summary of the number of tremors in terms of the depth of occurrence.

Figure 7. Interferograms (before unwrapping) calculated for four seismic events. Calculation results for paths (left) ascending (No. 73) and (right) descending (No. 22). The tectonic faults are marked on the maps with dashed lines: the Rudna Glowna fault (NW-SE) and the Biedrzychowa fault (SW-NE).

Figure 8. The modeling result for the analyzed mining tremors. Both vertical and horizontal displacements are shown for each event.

Figure 9. Line-of-Sight (LOS) displacements calculated based on the SBAS method. Successive curves correspond to 6 day increments of LOS displacements. In all cases, the calculated profiles cross the subsidence troughs in which permanent displacements after mining tremors from 2016, 2017 and 2019 were observed. The areas marked in red represent the size of the LOS displacement increases observed after the mining tremors occurred. The presented SBAS analysis covered a period from 15 November 2014 to 26 October 2020.

Figure 10. Comparison of the displacement values obtained before (black line) and after (green line) applying the tropospheric correction for the LGCB area. Maps show vertical displacements before and after correction, as well as E-W displacements before and after correction. Histograms show LOS displacements for paths (ascending and descending) before and after tropospheric correction.

Figure 11. Profiles (E-W cross-section) of the deformation indicators for modeled (dashed lines) and observed (solid lines) displacements caused by the analyzed mining tremors. For each analyzed tremor the upper profile represents subsidence (red) and horizontal shift (blue), while the lower profile shows slope (red) and strain (blue).

Table 1

Comparison of some induced seismic tremors around the world with regards to the type of exploitation. The table also includes examples of how synthetic aperture radar (SAR) data are used to facilitate the detection of ground surface changes due to induced seismic tremors.

Event Location	Date of Event	MaximumMagnitude [M${}_{\mathit{W}}$]	Type of Event	Examples in Literature	Examples of SAR Data Usage
Pawnee,Oklahoma, USA	3 September 2016	5.8	Water injection	[12]	[13,14,15]
Raton Basin,Colorado andNew Mexico	23 August 2011	5.3	Water injection	[16,17]	-
Bachatsky,Kuzbass, Russia	18 June 2013	6.1	Mining (coal)	[18]	-
Rudna,Poland	15 September 2018	4.8	Mining (copper)	-	[19]
Saar (Primsmulde),Saarland, Germany	23 February 2008	4.0	Mining (coal)	[20,21]	-
Pohang (PX-2),South Korea	15 November 2017	5.5	Geothermal	[22,23]	-
Fashing,Texas, USA	20 October 2011	4.8	Gas extraction	[24]	-
Selemo andLesedi pilot pods,Botswana	3 April 2017	6.3	Coal Bed Methane	[25]	[25]

Table 2

Indicators characterizing subsidence trough for points with known coordinates [57].

Indicator		Equation	Unit	Remarks
	Vertical	$w_i=H_i^a-H_i^p$	mm	Subsidence
Displacements	Horizontal	$u_i=\sqrt{\Delta X_{i_{a,p}}^2+\Delta Y_{i_{a,p}}^2}$	mm	Shift
Tilt		$T_{i,i-1}=\frac{w_i-w_{i-1}}{l_{i,i-1}}$	mm/m	Subsidence effect
Curvature		$K_{i-1,i+1}=\frac{T_{i,i+1}-T_{i-1,i}}{l_r}$	km${}^{-1}$	Tilt effect
Strain		${\epsilon }_{i,i-1}=\frac{l_{i,i+1}^a-l_{i,i+1}^p}{l_{i,i+1}^p}$	mm/m	Shift effect

a is the result of the current measurement, p is the result of the previous measurement (periodic change) or the initial measurement (total change), i − 1, i; and i + 1 mean three points adjacent to each other, H_i is the height of the i-th measurement point, l_i,i+1 is the distance between the measured neighboring points, and l_r is the average length of two adjacent sections of the observation line $l_r=\frac{l_{i,i-1}+l_{i,i+1}}{2}$.

Table 3

Scene-acquisition data: basic information.

		Ascending, Path no. 73		Descending, Path no. 22
Date and Time of Event	Strength [ M${}_{\mathbf{W}}$ ]	Master Date and Time	Slave Date and Time	Master Date and Time	Slave Date and Time
29 January 2019 12:53:44.3 a.m.	3.7	7 January 2019 04:43:36 p.m.	10 February 2019 04:43:48 p.m.	26 January 2019 05:09:02 a.m.	1 February 2019 05:08:27 a.m.
26 December 2017 11:15:30.2 a.m.	3.6	23 December 2017 04:42:51 p.m.	29 December 2017 04:43:33 p.m.	26 December 2017 05:08:57 a.m.	1 January 2018 05:08:21 a.m.
Events location: LGCB region, Poland; type: induced tremors
7 December 2017 05:42:50.0 a.m.	3.3	23 December 2017 04:42:51 p.m.	29 December 2017 04:43:33 p.m.	2 December 2017 05:08:58 a.m.	8 December 2017 05:08:22 a.m.
29 November 2016 08:09:41.1 p.m.	3.4	28 November 2016 04:43:19 p.m.	10 December 2016 04:43:19 p.m.	19 November 2016 05:08:02 a.m.	1 December 2016 05:08:01 a.m.

Dates and times are described in the following format: dd/mm/yyyy and hh:mm:ss (UTC). LGCB: Legnica-Glogow Copper Belt.

Table 4

Source parameters (moment tensor, source components and focal mechanism) of the selected seismic events.

	Coordinates			MT Decomposition			Mechanism
Date and Time of Event	Lat	Long	Depth	V(Explosive)	DC	CLVD(P-axis)	Plot
	(deg)	(deg)	(km)	(%)	(%)	(%)
29 January 2019 12:53:44 a.m.	51.5110	16.1197	0.8	32.7	8.8	58.5	[Image omitted. Please see PDF.]
26 December 2017 11:15:30 a.m.	51.5088	16.1065	0.7	30.3	9.1	60.5	[Image omitted. Please see PDF.]
7 December 2017 05:42:50 a.m.	51.5008	16.1021	0.9	26.6	16.4	56.9	[Image omitted. Please see PDF.]
29 November 2016 08:09:41 p.m.	51.5145	16.1573	0.8	26.0	21.7	52.2	[Image omitted. Please see PDF.]

Dates and times are described in the following format: dd/mm/yyyy and hh:mm:ss (UTC). CLVD: compensated linear vector dipole; DC: double couple forces.

Table 5

Source parameters determined from DInSAR results.

		Yang Models			Mogi Model
		29 November 2016	7 December 2017	29 January 2019		26 December 2017
Parameters	Depth, [km]	0.8	0.9	0.8	Depth, (km)	0.8
	Major axis, (km)	0.24	0.19	0.33	Radius (km)	0.06
	Minor axis, (km)	0.08	0.07	0.04	$\Delta P$$\times {10}^6$, (Pa)	−672
		$\Delta P$$\times {10}^6$, (Pa)	−282	−521	−709
	Strike, (deg)	33	313	195
	Plunge, (deg)	18	21	7