Tracking Deformation Processes at the Legnica Glogow Copper District (Poland) by Satellite InSAR—I: Room and Pillar Mine District

: Mining exploitation leads to slow or rapid ground subsidence resulting from deformation until the collapse of underground post-mining voids following excavation activities. Satellite SAR interferometry capabilities for the evaluation of ground movements allows the monitoring of intensive surface mine subsidence and can provide new knowledge about the risks in the mining industry. This work integrates both conventional and advanced Differential SAR Interferometry (DInSAR) to study the ground subsidence in the Legnica Glogow Copper District (LGCD, Poland) by processing about 400 Sentinel-1 images from October 2014 to April 2019. Even without ﬁeld data and information on past and ongoing excavation activities, the DInSAR approach allowed us to identify 30 troughs of subsidence, ranging from 500 m to 2.5 km in diameter, which in some cases, took place several times during the analyzed time span. The cumulative subsidence in 4 years and 7 months exceeds 70 cm in several zones of the LGCD. The sub-centimetric precision achieved by advanced analysis (A-DInSAR), allowed us to monitor the real extent of the mining inﬂuence area on the surface, with deformation velocities of up to 50 mm/year. The ground deformation detected at LGCD can be due to both mining-induced tremors and roof subsidence above the underground excavation rooms. As deformations do not occur concurrently with tremors, this can be related to excavation activities or to degradation of abandoned mines.


Introduction
As an active mining district, the Legnica Glogow Copper District (LGCD, Poland) represents a dynamic environment that evolves very quickly in time, following the ore exploitation activities in the room-and-pillar mines. The Lower Silesian region is widely exploited for copper extraction, and the LGCD mines represent an important source of ore managed by the KGHM (Polish acronym for Copper Mine and Mill Company), a mining company owned by the Polish government. The whole district is a risk-prone area that can be affected by subsidence and stability issues of the underground chambers and excavation fronts. Therefore, the active mining districts need accurate monitoring for preserving the safety of operators during the works and of nearby urban settlements.
The monitoring of mining subsidence by traditional techniques, such as field surveys, leveling, total stations, or GNSS, is labor-intensive and time-consuming [1,2] and it is usually limited to areas of a few square kilometers. Traditional techniques are indeed conducted on a point-by-point basis with limited spatial coverage, while the ground subsidence generally affects large areas [3]. Satellite remote sensing technologies such as Synthetic Aperture Radar (SAR) Interferometry (InSAR), can overcome this limitation, allowing us to analyze wide areas and obtain information from the past by means of data archives. InSAR has gained more and more utilization in mines monitoring, providing ground deformation measurements with sub centimeter accuracy [4][5][6][7][8]. This technique is very effective for structural health monitoring purposes involving structures like dams (e.g., [9][10][11] and Cultural Heritage degradation (e.g., [12])). Specifically, Advanced Differential SAR Interferometry (A-DInSAR) can provide information about the ongoing and past deformations with millimeter accuracy and can be used for many different applications in the field of engineering geology, e.g., subsidence and settlements [13,14], landslides (e.g., [15,16]), glaciers, and dynamics of permafrost [17,18].
In underground mines of copper ores, the most common and dangerous threats are represented by seismic tremors and rock bursts. In the room-and-pillar system, the deposit bed is cut (using blasting technique) into rooms and belts, which create technological pillars. The operating pressure causes damages to these pillars [19], which can evolve following a rock burst. An intact rock mass can therefore be affected by seismic tremors induced by mining work in an area that was previously aseismic. Tremors are the result of energy release due to the rock mass overstressing process in the pillars (the rock burst phenomenon). The related fracturing causes subsidence that affects not only excavations and underground infrastructure, but also the terrain surface. Therefore, the superficial ground deformation phenomena that can occur at active room and pillar mining districts can be related to both (i) land subsidence above the underground excavation rooms, produced at the advancing working face in the excavation direction (e.g., [20]) or due to yielding in the abandoned mines (e.g., [20]), and (ii) mining-induced tremors, related to post-pillar-failure disasters [19,20]. The effect of mining-induced tremors also results in a subsidence phenomenon, but it could take place in peculiar ways and times (e.g., [21]).
The (i) land subsidence due to underground ore exploitation is well known in literature and it can range from small-scale collapses to regional settlements. Its mechanism depends on the geological characteristics, the mechanical properties of the rocks, and the applied extraction technology [1,3]. During the excavation activities, in general the subsidence basin is expanding along the advancing direction of the working face [22,23] and the deformation rate increases consistently with the advancing distance, until it stabilizes at a critical value. The support provided by the reserved pillars decreases the deformation [22]. Moreover, especially in abandoned areas, the collapse of a single pillar affected by the coupled effects of weathering and stress can results in ground subsidence at the surface. Cascading pillar settlements may occur if the neighboring pillars cannot carry the additional load, or if the resulting roof span is too large for effective load transfer [20], and may occur in a couple of hours or decades after mine abandonment. This makes it difficult to establish the exact cause-effect relationship and the timing of the ground subsidence above the mines by geodetic measurements like InSAR.
Regarding mining-induced tremors (ii), an intense seismicity has been documented at LGCD (e.g., [24,25]): according to the local mining monitoring system and the national seismological survey carried out by the Institute of Geophysics in the Polish Academy of Sciences (IGF PAS), more than 2000 events with magnitudes between 0.9 and 3.5 are recorded yearly in this district [25]. The mining-induced earthquakes can cause the formation of large subsidence troughs occurring rapidly as a result of the release of seismic energy, with dimensions of a few square kilometers and total subsidence of the order of about 10 centimeters [26].
A-DInSAR algorithms, such as Persistent Scatterers Interferometry (PSI), can provide precise results regarding the relatively slow displacement, but this multitemporal technique has drawbacks in measuring fast deformations [27], due to ambiguity in the unwrapping procedure. Ng et al. [28] suggested the consecutive DInSAR method with a shorter temporal baseline as the most appropriate tool for fast terrain motions detection. A solution to Land 2021, 10, 653 3 of 20 deal with a strong subsidence in the center of the subsidence trough is the combination of two techniques: conventional (DInSAR) for the fast motions and advanced (A-DInSAR) for the residual subsidence [29]. Because traditional differential interferometry exploits image pairs, the accuracy of this technique is limited by factors related to spatial and temporal decorrelation, signal delays due to the atmospheric artifacts, and orbital or topographic errors [30]. Despite these problems, many scientists have successfully applied the consecutive DInSAR approach to the monitoring of fast-mining-related subsidence [3][4][5][6][7][8]29,31].
This paper focus on the deformation phenomena that affect the room-and-pillar LGCD copper mine ( Figure 1) using purely satellite SAR interferometry, and in particular combining both Advanced Differential SAR Interferometry (A-DInSAR) and Differential SAR Interferometry (DInSAR). Even without detailed information on current excavation activities, and without the support of field data, this work highlights the capabilities of SAR satellite monitoring of ground subsidence resulting from mining. Such a deformation phenomenon could introduce significant geological risk in the densely populated area (Lower Silesian voivodeship, south-western Poland) where several urban centers and critical infrastructure like theŻelazny Most Tailings Dam are located. The A-DInSAR analysis of this important storage facility, adjacent to the LGCD area, is described in a dedicated paper (Mazzanti et al., accepted).
Land 2021, 10, x FOR PEER REVIEW 3 of 22 combination of two techniques: conventional (DInSAR) for the fast motions and advanced (A-DInSAR) for the residual subsidence [29]. Because traditional differential interferometry exploits image pairs, the accuracy of this technique is limited by factors related to spatial and temporal decorrelation, signal delays due to the atmospheric artifacts, and orbital or topographic errors [30]. Despite these problems, many scientists have successfully applied the consecutive DInSAR approach to the monitoring of fast-mining-related subsidence [3][4][5][6][7][8]29,31]. This paper focus on the deformation phenomena that affect the room-and-pillar LGCD copper mine ( Figure 1) using purely satellite SAR interferometry, and in particular combining both Advanced Differential SAR Interferometry (A-DInSAR) and Differential SAR Interferometry (DInSAR). Even without detailed information on current excavation activities, and without the support of field data, this work highlights the capabilities of SAR satellite monitoring of ground subsidence resulting from mining. Such a deformation phenomenon could introduce significant geological risk in the densely populated area (Lower Silesian voivodeship, south-western Poland) where several urban centers and critical infrastructure like the Żelazny Most Tailings Dam are located. The A-DInSAR analysis of this important storage facility, adjacent to the LGCD area, is described in a dedicated paper (Mazzanti et al., accepted).

Geological Setting
The Legnica-Głogow Copper District (LGCD) is located in the Lower Silesia region (SW Poland), at the southern termination of the Fore-Sudetic Monocline. The latter is composed of metamorphic rocks and volcanogenic-sedimentary beds. It is represented by a thick epi-Hercynian sequence of sediments and rocks of the Lower and Upper Permian, which were respectively dated to the Rotliegende and Zechstein ages (Figure 2, [21]). The tectonic units that occur in the mining area are represented by: (i) the Proterozoic-Palaeozoic metamorphic basement, (ii) the Permo-Triassic complex of sedimentary

Geological Setting
The Legnica-Głogow Copper District (LGCD) is located in the Lower Silesia region (SW Poland), at the southern termination of the Fore-Sudetic Monocline. The latter is composed of metamorphic rocks and volcanogenic-sedimentary beds. It is represented by a thick epi-Hercynian sequence of sediments and rocks of the Lower and Upper Permian, which were respectively dated to the Rotliegende and Zechstein ages (Figure 2, [21]). The tectonic units that occur in the mining area are represented by: (i) the Proterozoic-Palaeozoic metamorphic basement, (ii) the Permo-Triassic complex of sedimentary rocks, which is smoothly dipping to the NE, and above the latter, and (iii) a more or less horizontal Cenozoic sediments sequence [32]. The copper ore deposits belong to the Permo-Triassic complex of sedimentary rocks and in particular to the Kupferschiefer Formation and to the basal part of the Zechstein Limestone, and the mines are localized at the boundary between the platform and the basin facies of the Zechstein Limestone ( [33]; Figure 2a rocks, which is smoothly dipping to the NE, and above the latter, and (iii) a more or less horizontal Cenozoic sediments sequence [32]. The copper ore deposits belong to the Permo-Triassic complex of sedimentary rocks and in particular to the Kupferschiefer Formation and to the basal part of the Zechstein Limestone, and the mines are localized at the boundary between the platform and the basin facies of the Zechstein Limestone ( [33]; Figure 2a,b). Stratiform mineralization founded within sandstones, shales, and dolomites, occurs in different amounts from mine to mine and is located at depths between 600 m and 1200 m (www.mining-technology.com, accessed on 15 May 2021). The roof of the LGCD mines consists of sub-horizontal limestone or dolomite strata 200 m thick, while anhydrite and gypsum layers, parallel to the bedding, lie above ( Figure 2b). The ore is located above the basement rocks at the intersection with two NW-SE and NNE-SSW striking fault systems that control the regional distribution of the metal and the content of mineralization and its geometry [34]; the ore is tectonically disturbed by a dense grid of fractures and faults [35,36]. The dominant fault system in the mining area coincides with the fault zone of Rudna and Lubin (Figure 2a). These two faults represent parallel structures striking Stratiform mineralization founded within sandstones, shales, and dolomites, occurs in different amounts from mine to mine and is located at depths between 600 m and 1200 m (www.mining-technology.com, accessed on 15 May 2021). The roof of the LGCD mines consists of sub-horizontal limestone or dolomite strata 200 m thick, while anhydrite and gypsum layers, parallel to the bedding, lie above ( Figure 2b). The ore is located above the basement rocks at the intersection with two NW-SE and NNE-SSW striking fault systems that control the regional distribution of the metal and the content of mineralization and its geometry [34]; the ore is tectonically disturbed by a dense grid of fractures and faults [35,36]. large stratigraphic throws, and the Rudna fault is considered the most important due to the potential seismic hazard ( Figure 2).

Materials and Methods
InSAR (Interferometric Synthetic Aperture Radar, [37]) is based on the comparison between RADAR images acquired in different times over the same area. Differential SAR Interferometry (DInSAR) has been successfully applied in the investigation of different ground deformation processes and is mainly characterized by impulsive or intermittent behavior. Indeed, DInSAR provides good results for the investigation of displacement characterized by wide, spatially smooth deformation processes, such as coseismic and postseismic deformations [38][39][40], volcanic deformation processes [41,42], or ice and glacier dynamics [17,43] or landslide detection [44][45][46]. Furthermore, mining-related deformations phenomena have also been detected and measured by DInSAR (e.g., [3,29]).
On the other hand, the Advanced Differential SAR Interferometry (A-DInSAR) approach consists of the combination of information from a large number of SAR images, allowing us to derive the temporal evolution of displacements of objects on the ground during the investigated period. Successful applications to settlement/subsidence-related issues have been carried-out through an A-DInSAR approach in the last few years [3,13,14,29,47,48].
In this study, Sentinel-1 SAR images were analyzed using a hybrid approach based on the standard Differential Interferometry (DInSAR) and Multi-Image techniques.
As part of this process, 411 C-band and Sentinel-1 SAR scenes (in both ascending and descending orbital geometries) were acquired from the ESA archives for a time span of about 5 years from October 2014 to April 2019. We used Interferometric Wide beam mode in Single Polarization (IW) SAR images, characterized by a spatial resolution of approximately 5 × 20 m, a wavelength of 5.6 cm (C-band), and a 250 km wide images (Table 1).

Multi-Temporal DInSAR Analysis
The DInSAR technique is based on the analysis of the phase variations by coupling 2 SAR images collected at different times (namely "interferograms"). The differential interferograms were computed to consider both temporal and spatial coherence. For this reason, we performed a multi-temporal DInSAR analysis based on Minimum Spanning Tree graph (MST, [49]), creating a framework for optimal interferometric pairs selection. MST involve a set of interferograms from SAR stack corresponding to the minimum best coherent graph, which were obtained by connecting all the images of the dataset. A total amount of 411 interferograms between couples of images in the same orbital geometry spanning for the entire monitoring period were processed. The availability of images both in ascending and descending orbital geometry allowed us to have a redundant and independent source of information.
A key step to achieve reliable results about displacement occurrences was the detection and subtraction of the topographic component on differential interferograms. The Shuttle Radar Topography Mission (SRTM) 1 arc-sec (30 m resolution) DEM, provided by the US Geological Survey, was used to remove the topographic component from the generated interferograms.
The generated interferograms were unwrapped [50,51] to create deformation maps which show displacement values continuously in space. The Goldstein filter [52] was applied to support the Minimum Cost Flow (MCF) phase unwrapping procedure. The where λ is the wavelength used by the C-band sensor. The displacement measured by the unwrapped interferograms was cumulated over time to obtain a cumulative displacement map, with centimetric accuracy. In particular, the cumulative surface displacement was computed using the contribution of interferograms generated from consecutive pairs of SAR images.

Advanced DInSAR Analyses (A-DInSAR)
Advanced Differential Interferometric SAR (A-DInSAR) analysis was performed to acquire information regarding displacements that affected the ground surface over the LGCD area since 2014.
A large-scale study was performed on an area of about 700 km 2 , carrying out a conventional PSI approach in accordance with Ferretti et al. [53,54]. Specifically, to reduce the decorrelation effects, all images were related to a single master image, which was selected by considering the perpendicular and temporal baselines. After all images were co-registered to the master image, the maps of radar parameters were generated: the reflectivity map (multi-temporal amplitude value for each pixel) and the amplitude stability index map (coefficient of the amplitude variation). These maps were used to estimate the quality of the persistent scatterer candidates (PSCs) as a first step of the PSI workflow. A network of PSCs was created to estimate preliminary height and displacement parameters. This step is needed to detect and remove the atmospheric phase screen (APS) starting from the residual phase components. For each persistent scatterer (PS), the velocity values along the Line of Sight (LOS), displacement time series, and height were computed. It is worth mentioning that (i) displacements were measured along the satellite LOS as the path between the satellite sensor and the target (e.g., look angle in Table 1), and that (ii) these values of displacements were relative to a reference point identified in a stable region. At the end of the PSI workflow, the more reliable PSs were selected by applying a high temporal coherence threshold (>0.6). For such PSs, the time series of displacements has been provided.
Lastly, the Up-Down and East-West displacement rates were retrieved by combining the ascending and descending orbital geometry [58]. An application for vector decomposition owned by NHAZCA S.r.l. was used for this purpose, according to the following steps: (1) the study area was discretized in a regular grid of hexagonal cells; (2) when the cells contained at least one PS both from the ascending and the descending geometry, the vertical and horizontal velocity rates could be achieved; (3) the ascending and descending LOS velocities were projected along the vertical and horizontal directions (up-down and east-west) by a vector decomposition. The horizontal component along the north-south direction cannot be measured because the satellites fly in a near-polar orbit approximately perpendicular to the LOS.  Figure 4. The time series were retrieved averaging the time series for 3 neighboring PSs, so that the observed trend is more representative for the areas, without any isolated outliers.

Results
Land 2021, 10, x FOR PEER REVIEW 7 of 22 been observed in correspondence with Tarnowek, near the tailing basin and east of Polkowice. Regarding the horizontal deformations (east-west direction), the maximum velocity values were up to 46 mm/year toward east (blue colors in Figure 3a), observed in the city of Polkowice, and −28 mm/year toward west (orange colors in Figure 3a). The east-west and up-down time series of displacement for the Polkowice, Tarnowek, and Moskorzyn areas are shown in Figure 4. The time series were retrieved averaging the time series for 3 neighboring PSs, so that the observed trend is more representative for the areas, without any isolated outliers.  been observed in correspondence with Tarnowek, near the tailing basin and east of Polkowice. Regarding the horizontal deformations (east-west direction), the maximum velocity values were up to 46 mm/year toward east (blue colors in Figure 3a), observed in the city of Polkowice, and −28 mm/year toward west (orange colors in Figure 3a). The east-west and up-down time series of displacement for the Polkowice, Tarnowek, and Moskorzyn areas are shown in Figure 4. The time series were retrieved averaging the time series for 3 neighboring PSs, so that the observed trend is more representative for the areas, without any isolated outliers.  deformation of each zone is gradual, i.e., highly deforming areas are surrounded by zones that show minor and gradual deformation.
To investigate the deformation pattern observed in the study area in detail, the DInSAR (Differential SAR Interferometry) approach was used, with the analysis and interpretation of 400 interferograms relating to the time period from October 2014 to April 2019. In the interferograms, the areas affected by deformation (with rates in the order of some centimeters for month) are represented by the so-called interferometric fringes, i.e., bands of color from blue to red, which corresponds to a ground displacement of 2.8 cm each (for C-band sensors). The deformation areas are characterized by concentric fringes, with sub-circular or ellipsoidal shapes, with dimensions ranging from 500 m to 2.5 km in diameter. Although the observed fringes represent a LOS displacement (moving away from the sensor), the latter was considered here as vertical deformation, namely subsidence, neglecting the possible occurrence of minor horizontal components. Figure 5 shows an interferogram with temporal baseline of 48 days that represents several subsidence troughs occurring sparsely throughout the mining area. The zones affected by intense deformation were mapped in all the available interferograms for a total of 30 subsidence troughs in the whole analyzed time span. The velocity maps obtained by the A-DInSAR approach (Figure 3a,b) highlight the presence of some widespread deformations above the room and pillar mines of the LGCD, showing a complex and discontinuous pattern in space. Furthermore, the trend of the deformation of each zone is gradual, i.e., highly deforming areas are surrounded by zones that show minor and gradual deformation.
To investigate the deformation pattern observed in the study area in detail, the DInSAR (Differential SAR Interferometry) approach was used, with the analysis and interpretation of 400 interferograms relating to the time period from October 2014 to April 2019. In the interferograms, the areas affected by deformation (with rates in the order of some centimeters for month) are represented by the so-called interferometric fringes, i.e., bands of color from blue to red, which corresponds to a ground displacement of 2.8 cm each (for C-band sensors). The deformation areas are characterized by concentric fringes, with sub-circular or ellipsoidal shapes, with dimensions ranging from 500 m to 2.5 km in diameter. Although the observed fringes represent a LOS displacement (moving away from the sensor), the latter was considered here as vertical deformation, namely subsidence, neglecting the possible occurrence of minor horizontal components. Figure 5 shows an interferogram with temporal baseline of 48 days that represents several subsidence troughs occurring sparsely throughout the mining area. The zones affected by intense deformation were mapped in all the available interferograms for a total of 30 subsidence troughs in the whole analyzed time span.  it is possible to count fringes and to estimate several centimeters of subsidence. It was observed that up to several fringes can be counted before the interferogram is decorrelated in the central part, probably due to the magnitude of the displacement beyond the detection threshold of the sensor. The fringes at some subsidence troughs occur on several subsequently analyzed interferograms, suggesting a continuous subsidence over time. For example, in both of the 48 days interferograms of Figure 6a,b, some of the subsidence troughs show similar deformation fringes during different time periods. Figure 5. Sentinel 1 C-band differential interferogram, dates: 2015/01/14-2015/03/03; white lines represent LGCD mine district borders. The perimeter of the subsidence troughs was mapped during the whole analyzed time span (October 2014-April 2019). Figure 6 represents seven interferograms of a subset of the study area (rectangle in Figure 5) located west of Żelazny Most dam, near Polkowice town. In the interferograms, it is possible to count fringes and to estimate several centimeters of subsidence. It was observed that up to several fringes can be counted before the interferogram is decorrelated in the central part, probably due to the magnitude of the displacement beyond the detection threshold of the sensor. The fringes at some subsidence troughs occur on several subsequently analyzed interferograms, suggesting a continuous subsidence over time. For example, in both of the 48 days interferograms of Figure 6a,b, some of the subsidence troughs show similar deformation fringes during different time periods. On the other hand, some deformation zones show a discontinuous pattern over time, and in several cases, they disappear suddenly in subsequent interferograms or do not maintain the same intensity of subsidence, showing relatively intermittent or impulsive behavior. For example, in Figure 6g, a subsidence trough shows four fringes in the interferogram that spans from 26 January 2019 to 2 February 2019, which roughly correspond to about 11.2 cm of displacement in 6 days, but fringes were totally lacking in several previous interferograms (Figure 6c,e,f).
The intensity of the subsidence could be quantitatively estimated by the analysis of unwrapped interferograms that represent deformation maps with centimetric accuracy and are effective for analyzing the displacement, as was recently demonstrated by Borg et al. [21], Dai et al. [49], and Sui et al. [59]. To quantitatively estimate the total displacement that occurred during the whole analyzed time interval of about 4 years and 7 months (October 2014-April 2019), the unwrapped interferograms (in ascending geometry) were cumulated in a map showing the cumulative displacement (Figure 7). Such a map shows that the deformation persists over time in the same areas (yellow to red colors in Figure 7), which coincide with the subsidence troughs. For 10 different zones, the maximum cumulative subsidence exceeds 70 cm in about 4 years and 7 months. Where some subsidence troughs are very close to each other, the surface affected by subsidence includes quite large areas, which were mapped in Figures 7 and 8. These zones are localized mainly in the northern and central parts of LGCD, at Sieroszowice and Rudna mines, respectively.      Figure 8a) that represent the cumulative displacement profiles above the rock excavations. For each selected subsidence trough, the cross sections were traced in an orthogonal manner, and their variable lengths include the area affected by surface deformation.
The profiles allowed us to appreciate the intensity and the spatial pattern of the cumulated deformation during the analyzed time span. The profiles of subsidence trough A show a slightly asymmetrical shape, while the displacement curves in C have a symmetrical shape. On the other hand, the displacement curve of subsidence trough B in NS direction shows that subsidence pattern continues towards the south due to a larger area that is affected by moderate ground deformation.

Discussion
The deformation pattern shown by the A-DInSAR velocity maps (Figure 3a,b) indicates several localized deformation zones above the room and pillars mines of the LGCD. However, large areas are devoid of measurement points (PSs) despite including urbanized zones, where many PSs are normally detected by the PSI technique. Conversely, zones showing intense deformation, with a high density of PSs, surround these areas where PSs are lacking. Such a particular PSs distribution can be due to the occurrence of too rapid displacement in the areas devoid of PSs that cannot be properly investigated by the A-DInSAR technique, as was argued by some authors (e.g., [5,29]). Ambiguity in the unwrapping procedure results in the limitation of the maximum differential displacement rate between two neighboring points that can be measured by PSI. This value is limited to λ/4 over the revisiting time between two consecutive SAR images, where λ is the wavelength of the radar signal (λ = 5.6 cm for C-band satellites). With a revisit time of six days (minimum temporal baseline of the interferograms), the theoretical maximum displacement rate detectable by PSI avoiding ambiguity in the unwrapping is 85 cm/year. This value is reached only if (i) all SAR images every 6 days are available and (ii) if the observed deformation is constant and equal to 1.4 cm every 6 days. Higher values would result in decorrelation and would be lost in the total count. The theoretical maximum displacement rate detectable by PSI decreases when the temporal baseline of the interferogram increases, or when the deformation process is not constant. These are the reasons why it has not been possible to measure the displacement in the center of the subsidence through using PSs, despite the fact that the cumulative displacement (Figure 7) is lower than this theoretical maximum rate.
Despite this limitation, the sub-centimetric precision achieved by PSI analysis, allowed us to monitor the residual subsidence and to map the more moderate ground displacements surrounding several intense subsidence areas of limited sizes detected by the DInSAR classical approach. The latter therefore allowed us to identify the subsidence troughs that were previously recognized in the literature in LGCD and in other rooms and pillar mines [26,60]. In agreement with some authors (e.g., [28]), the results confirmed that the analysis of interferograms (Figures 5 and 6) represents an appropriate tool for detecting intense and fast terrain motions where the advanced multitemporal technique is no longer efficient.
In Figure 9, the complementarity between the two different SAR interferometry techniques (A-DInSAR and DInSAR) is highlighted: colored points represent the PSs in an up-down direction and the black and white image represents the cumulative DInSAR displacement map of Figure 7 in shades of gray. In Figure 9c, it is worth noting that in the areas not covered by PSs, several sub-circular or ellipsoidal fringes occur. This is clear in Figure 9b where a subset of interferogram near Polkowice is shown: where the PSs are lacking, the interferogram shows ellipsoidal fringes, which correspond to the areas with the highest displacement rate. This residual subsidence located around the areas affected by intense deformation is well mapped by PSs distribution and decreases gradually by up to 1.5 mm/year (considered as a stability threshold) about 1.5 km away from the center of the subsidence troughs. Land 2021, 10, x FOR PEER REVIEW 14 of 22 In the LGCD region, 30 different ellipsoidal fringes covering relatively small areas indicate troughs of subsidence, which in some cases took place several times in the same area during the analyzed time span. The location of these active subsidence areas within the perimeter of the mines suggests a direct connection between the mining activity and the superficial deformation phenomena (Figures 5 and 7). This hypothesis is supported by a previous work of Popiolek et al. [61] that carried out a pioneering application of the InSAR technique for the monitoring of LGCD areas using a single interferogram processed from two ERS-1 images of January and March 1994. These authors argued that the origin of the 26 detected subsidence troughs was directly connected with the mine exploitation by comparing the analysis of the interferogram and the mining exploitation development for the analyzed time.
The behavior over time of the subsidence phenomena is rather complex to analyze: 13 troughs continue to undergo subsidence more or less constantly, while 17 areas are affected by deformation pulses that take place for short time intervals of the order of a few days, and that can be repeated several times during the analyzed years ( Figures 5  and 10). In the LGCD region, 30 different ellipsoidal fringes covering relatively small areas indicate troughs of subsidence, which in some cases took place several times in the same area during the analyzed time span. The location of these active subsidence areas within the perimeter of the mines suggests a direct connection between the mining activity and the superficial deformation phenomena (Figures 5 and 7). This hypothesis is supported by a previous work of Popiolek et al. [61] that carried out a pioneering application of the InSAR technique for the monitoring of LGCD areas using a single interferogram processed from two ERS-1 images of January and March 1994. These authors argued that the origin of the 26 detected subsidence troughs was directly connected with the mine exploitation by comparing the analysis of the interferogram and the mining exploitation development for the analyzed time.
The behavior over time of the subsidence phenomena is rather complex to analyze: 13 troughs continue to undergo subsidence more or less constantly, while 17 areas are affected by deformation pulses that take place for short time intervals of the order of a few days, and that can be repeated several times during the analyzed years ( Figures 5 and 10). In the LGCD, in which the room and pillar method with roof deflection and maintaining the central part of the mining field was used [62], the assessment of subsidence timing and impact was difficult to define, particularly when specific information on the excavation activity is not available, as in the case of this work. A criterion to assess the subsidence impact was based on empirical observation of ground movements and damage to the buildings above or near the mines [63].
To interpret the origin of the deformation phenomena that occur at the LGCD active room and pillar mining district, two possible causes were considered in this study: (i) subsidence above the underground excavation rooms related to both the advancing working face in the excavation direction (e.g., [20]), or pillar failures in abandoned areas which are usually accompanied by sudden roof subsidence (e.g., [20]) or (ii) mining-induced tremors ("coseismic" subsidence e.g., [21]). The mining-induced seismic events in this area show similar characteristics of weak earthquakes, but with a lower depth of the hypocenters (up to 1 km). As mining-induced earthquakes are frequent at LGCD (e.g., [24]), a possible direct relationship between these tremors and the subsidence troughs detected by DInSAR was investigated, relying on the existing literature and analyzing the obtained results in detail.
The mining-induced seismicity and related rock burst hazard may be the function of several parameters: in LGCD, the principal factors influencing the seismic hazard level In the LGCD, in which the room and pillar method with roof deflection and maintaining the central part of the mining field was used [62], the assessment of subsidence timing and impact was difficult to define, particularly when specific information on the excavation activity is not available, as in the case of this work. A criterion to assess the subsidence impact was based on empirical observation of ground movements and damage to the buildings above or near the mines [63].
To interpret the origin of the deformation phenomena that occur at the LGCD active room and pillar mining district, two possible causes were considered in this study: (i) subsidence above the underground excavation rooms related to both the advancing working face in the excavation direction (e.g., [20]), or pillar failures in abandoned areas which are usually accompanied by sudden roof subsidence (e.g., [20]) or (ii) mining-induced tremors ("coseismic" subsidence e.g., [21]). The mining-induced seismic events in this area show similar characteristics of weak earthquakes, but with a lower depth of the hypocenters (up to 1 km). As mining-induced earthquakes are frequent at LGCD (e.g., [24]), a possible direct relationship between these tremors and the subsidence troughs detected by DInSAR was investigated, relying on the existing literature and analyzing the obtained results in detail.
The mining-induced seismicity and related rock burst hazard may be the function of several parameters: in LGCD, the principal factors influencing the seismic hazard level (measured by the number of registered seismic events and the total seismic energy released) are the regional geological setting, i.e., the effective mining thickness and the thickness of burst-prone roof rock strata, and the possible variation of tectonic stress intensity [36]. Limestones and anhydrites that lie above the ore store elastic strain energy that can be released through fracturing and rock bursts [19]. Moreover, the depth of the excavations (that varies from 600 m up to 1200 m in LGCD) and the retreat mining method of excavation can provide favorable conditions for energy accumulation [64].
In the literature, there are some studies that focus on the use of the InSAR technique for analyzing the superficial deformation caused by mining-induced tremors in the study area [8,21,26,65,66] and in other sites [65,67]. Some of the major induced seismic events that struck the mining panels of LGCD district were analyzed by using SAR images: Malinowska et al. [21] and Milczarek [64] analyzed the surface deformation processes caused by the earthquakes of 29 November 2016 by DInSAR and by the Small Baseline Subset technique (SBAS), respectively. Within about 2 days after the earthquake, the ground surface rapidly subsided, forming an asymmetrical-shape trough, with a maximum post-earthquake subsidence of about 80 mm. Surface deformations vanished until about the 7th day after the earthquake [21]. Hejmanowski et al. [26] analyzed the spatio-temporal distribution of ground movements caused by the two tremors that occurred in 2017 with magnitude of 4.7 and 4.8 using the InSAR technique. The tremor was preceded by few millimeters per day of surface movements caused by on-going mining activity. A total of 24 h after the tremor, the main subsidence of about 70-80 mm took place. In the subsequent 4-5 days, the surface movements stopped.
To examine a possible cause-effect relationship between tremors and subsidence pattern detected at LGCD, the space-time coincidence between earthquakes and subsidence events was investigated. As described by some authors (e.g., [26]), the onset of intense subsidence would occur immediately after the tremor, or after a maximum of one day. Figure 10 shows the distribution of epicenters of seismic events with M > 3 which occurred from October 2014 to April 2019, and the locations of the subsidence troughs detected by DInSAR that were classified according to their deformation behavior: impulsive or continuous over time. Although Figure 10 shows that several epicenters are located in correspondence with the subsidence troughs, in most cases the dates of tremors and subsidence events (fringes) did not match. The lack of time coincidence is highlighted in the example shown in Figure 11, which represents the location of a subsidence trough between Polkowice town and theŻelazny Most reservoir, where the cause-effect relationship between tremors and deformation pulses was not verified. The graph shows (i) the dates of 5 deformation pulses (quantified in LOS displacement in the histogram bars) which took place in this area and represented in 5 different interferograms and (ii) the dates of the major seismic events (black stars) that occurred within a radius of 2.5 km from this subsidence area. Deformations do not occur concurrently with earthquakes, and a proxy of the direct cause and effect relationship is lacking.
Furthermore, from the analysis of the cumulative displacement map (Figure 7), it was possible to argue that many areas reach a total subsidence of up to several tens of centimeters that can be achieved by a continuous subsidence over time, or by the sum of many subsidence pulses. Such kinds of deformation behavior cannot be easily described as direct consequences of tremors, but according to other authors, it can instead be related to excavation activities or to the abandonment of tunnels. For example, Przylucka et al. [5] used conventional DInSAR to measure a very fast mining-induced subsidence rate (up to 166 cm/year), showing that these deformations coincided with underground mining excavations. Other recent works [22,23,68] use conventional D-InSAR and other methods to retrieve a precise three-dimensional mining deformation field in the case of near-horizontal coal seam longwall mines in China. In this particular condition, they show that the surface movement basin has symmetrical morphological characteristics and expands along the advancing direction of the working face during excavation activities. Different geological conditions or the presence of large faults, folds, aquifers, and different topographic heights can strongly influence the surface movement, as does the mining and roof management methods [68]. Furthermore, from the analysis of the cumulative displacement map (Figure 7), it was possible to argue that many areas reach a total subsidence of up to several tens of centimeters that can be achieved by a continuous subsidence over time, or by the sum of many subsidence pulses. Such kinds of deformation behavior cannot be easily described as direct consequences of tremors, but according to other authors, it can instead be related to excavation activities or to the abandonment of tunnels. For example, Przylucka et al. [5] used conventional DInSAR to measure a very fast mining-induced subsidence rate (up to 166 cm/year), showing that these deformations coincided with underground mining excavations. Other recent works [22,23,68] use conventional D-InSAR and other methods to retrieve a precise three-dimensional mining deformation field in the case of near-horizontal coal seam longwall mines in China. In this particular condition, they show that the surface movement basin has symmetrical morphological characteristics and expands along the advancing direction of the working face during excavation activities. Different geological conditions or the presence of large faults, folds, aquifers, and different topographic heights can strongly influence the surface movement, as does the mining and roof management methods [68].
Moreover, sudden roof subsidence, pillar spalling, or floor heave can be caused by pillar failures, that can appear more frequently when width-height ratios of pillars are smaller than 3 for coal mines, less than 2 for non-metal mines, and less than 1 for metal mines [20,69]. Studies on the long-term stability of mineral pillars showed that weathering and stress will eventually result in pillar size reduction and pillar strength degradation over time [70,71]. The dissolution of minerals indeed may reduce the size of the pil- Moreover, sudden roof subsidence, pillar spalling, or floor heave can be caused by pillar failures, that can appear more frequently when width-height ratios of pillars are smaller than 3 for coal mines, less than 2 for non-metal mines, and less than 1 for metal mines [20,69]. Studies on the long-term stability of mineral pillars showed that weathering and stress will eventually result in pillar size reduction and pillar strength degradation over time [70,71]. The dissolution of minerals indeed may reduce the size of the pillars [72]. Because of the possible combination of weathering and stress effects, the constant pillar scaling rate is empirical and site-specific, and it is difficult to obtain reliable results on pillar stability or pillar life [20].

Conclusions
Both A-DInSAR and DInSAR analyses performed on C-band Sentinel 1A SAR images shed light on a complex pattern of ground subsidence at an LGCD mine district, demonstrating that accurate monitoring is very important for preserving the safety of nearby urban centers. Even without information on mining activities, nor field data, this work highlights the capabilities of SAR satellite monitoring of superficial deformation resulting from mining.
A-DInSAR was exploited using the PSI technique to detect a maximum displacement velocity up to 50 mm/year and providing a dense measurement points distribution. Given the typical limitations of the technique, the standard PSI approach failed to detect fast ground movements, which occur within the subsidence troughs. Despite this limitation, the sub-centimetric precision achieved by PSI analysis, allowed us to monitor the real extent of the mining influence area, i.e., the more moderate but still dangerous ground deformation surrounding the detected subsidence troughs. This residual subsidence is located around the areas affected by intense deformation and decreases gradually (below 1.5 mm/year, considered as a stability threshold) about 1.5 km away from the center of the subsidence troughs.
The conventional DInSAR approach allowed us to detect deformation fringes of the order of several centimeters per months, and to define the location, the shape, and rough velocity values of subsidence in the mining area. To quantitatively estimate the total displacement that occurred during the whole analyzed time interval of April 2014-October 2019, the unwrapped interferograms were cumulated to develop a cumulative displacement map. For 10 different zones, the maximum cumulative subsidence exceeded 70 cm in about 5 years.
Two possible causes were investigated to explain the ground deformation detected at LGCD: mining-induced tremors or roof subsidence above the underground excavation rooms. Most of the observed subsidence phenomena could not be considered as direct consequences of tremors, but rather could be related to sudden or progressive roof subsidence of the excavation rooms, which directly follows the excavation activity. This hypothesis is supported by Popiolek et al. [61], who highlighted correspondence between the InSAR results and the LGCD exploitation development for the analyzed time. Another cause could lie in the progressive deterioration of abandoned mines. Because of the combination of weathering and stress effects, the stability of a room and pillar mine area is empirical and site-specific, making it very difficult to make assumptions about the precise causes and timing of the deformation processes, and the lack of information about the excavation program performed at LGCD mines and the location of active mining fronts prevented us from analyzing each single subsidence troughs in detail.
The new generations of satellites ("SmallSats"; [73]) with very high-resolution sensors (up to <1 m) onboard and very short revisiting times (up to multiple times a day), are expected to expand the capabilities of the InSAR performances [74,75] and improve the monitoring of very fast deformation phenomena, which can affect a heavily anthropized environment that changes continuously and quickly, such as active mining areas. Funding: This research was funded by Sapienza Università di Roma, funding program: "Progetti di Ricerca Medi SAPIENZA 2018", grant number RM118164367B3FF7, title of the project: "Analisi dei rischi geologici interferenti con le grandi dighe attraverso l'uso integrato di dati telerilevati e realtà virtuale".