Mining and mineral extraction is useful for human development, social progression, and material prosperity. However, these activities strongly interact with surrounding environments and could produce meaningful negative effects, on the atmosphere, soil, and water (Bell et al. 2000). One of the main issues in many parts of the world is represented by the subsidence induced by underground mining (Villegas et al. 2011; Zhou et al. 2015; Przyłucka et al. 2015a, b; Vu Khac et al. 2021). Such effects are more harmful when the activities legal or not authorized are in correspondence or nearby urban areas (Bell et al. 2005; Prakash et al. 2010; Ji et al. 2011; Villegas et al. 2011; Marschalko et al. 2012). Subsidence can either coincide with mining activities or be delayed in response to the time-dependent deformation of rocks. Several studies determined the rate of expected subsidence employing physical or numerical modeling and analysis (Whittaker 1989; Singh et al. 1998; Alheib et al. 2001; Ye et al. 2016; Xie et al. 2020; Gazzola et al. 2021). In any case, to avoid possible environmental harm it is important to monitor areas subjected to underground mining activity. Mining-related subsidence is usually analyzed either by in-situ measurements, or by aerial photographs on specific sites (Armenakis 1983). As far as in-situ survey techniques concerns, they consist of: piezometers, useful to register pore pressure variations of the overburden strata (Guo et al. 2012); extensometers installed at different heights in boreholes located above the exploited panels: in this way, settlments are measured (Holla et al. 2000); and GNSS technique with benefits of accuracy and simultaneous 3-D positioning (Chen et al. 2014). However, these techniques are tedious and time-consuming, with small area coverage. Indeed, they have a cost of implementation that is significantly expensive for the management and the acquisition of information on large spatial and temporal scales of the entire mining area, thus providing localised and not distributed results.

Due to these limitations, remote sensing can be a cost-effective technique. Specifically, among the different types of remote sensing techniques, the Advanced Differential Interferometric Synthetic Aperture Radar (A-DInSAR—Franceschetti et al. 1992) proved to be a very powerful tool, being sensitive to sub-centimetric ground movements (Colesanti et al. 2006) and/or terrain displacements, induced by landslides (Wasowski and Bovenga 2022; Scifoni et al. 2016; Pappalardo et al. 2018; Giardina et al. 2019; Guerriero et al. 2019), earthquakes (De Novellis et al. 2018; Zhao et al. 2021), volcanoes (Foumelis et al. 2016; Casu et al. 2019) and valid to detect harm to structures and infrastructures (Milillo et al. 2018; Ullo et al. 2019; Pastor et al. 2019; Miano et al. 2021). In the last three decades, A-DinSAR has been applied also in mining areas to detect subsidence due to underground activity (Yue et al. 2011; Du et al. 2016; Ammirati et al. 2020; Pawluszek et al. 2020; Chen et al. 2021), monitoring the stability of tailings dam (Necsoiu et al. 2015; Gama et al. 2019; Ammirati et al. 2021), and to identify surface movements in open-pit mines (Paradella et al. 2015; Carlà et al. 2018). The scientific community has studied many cases of mining subsidence in coal districts (Salmi et al. 2017; Jing et al. 2018; Zingano et al. 2019). In these areas, A-DinSAR resulted to be an effective technique for ground movement surveys (Dong et al. 2013; Xu et al. 2020). In several studies, the subsidence maps generated from satellite images proved to be consistent with field observations (Ismaya et al. 2012; Samsonov et al. 2013; Pawluszek et al. 2020).

From this point of view, it is essential to not underestimate the potential effects of mining activities and constant monitoring the mineral resource exploitation and its effects on the ecological environment. Furthermore, mining monitoring is of greatest importance, not only for the scientific community but also for the mining companies, to reduce the potential effects and consequences.

In this work, an integrated approach based on geotechnical numerical modeling and A-DInSAR method was applied for investigating subsidence phenomena which occurred in a coal mining area affected by underground exploitation.

The innovativeness of this application lies in several aspects: (1) it is one of the first applications in the mining sector on Italian territory, (2) the time interval analyzed is about 10 years, thus allowing coverage of both the pre-, syn- and post-mining phases, (3) unlike previous works very high-resolution COSMO-SkyMed images were processed; in fact, there are few applications in the world that have used very high-resolution images for monitoring underground mining areas (Przyłucka et al. 2015a, b; Alam et al. 2022), finally, (4) it is one of the first combined applications, in the Italian mining sector, between a research institution and mining company with the objective of monitoring the rehabilitation of a dismissed mining area.

The study area is represented by a historic coal mining district located in southwest Sardinia, Italy, in the municipality of Nuraxi Figus (Carbonia-Iglesias province). The area occurs in the Carbosulcis S.p.A. mining concession called “Monte Sinni” (total extension of about 0.6 km²), where a subsidence phenomenon occurred between the years 2011 to 2014 (Tessitore et al. 2018). In the present study, we analyzed the remote sensing data, obtained in the framework of Not-Ordinary Plan of Environmental Remote Sensing (Piano Straordinario di Telerilevamento Ambientale—in Italian), funded by the Italian Ministry of Environment (Costantini et al. 2017; Di Martire et al. 2016), to correlate the 2011–2014 subsidence phenomenon with the underground mining activity and the excavation development, also using numerical modeling to reconstruct the kinematic evolution of ground displacements. Moreover, radar satellite data, derived from COSMO-SkyMed (CSK) datasets, were carried out by SUBSIDENCE software, which implements the Coherent Pixels Technique algorithm (Mora et al. 2003; Iglesias et al. 2015), for the period from 2013 to 2020, to understand the temporal evolution of ground surface deformation. The study allowed us to understand the cause-effect mechanism between extraction and ground deformation and confirmed that ground surface deformations can continue at lower rates a few years after the completion of mining activities.

In the Sulcis coal district, the mining activity started around 1850 through the exploitation of productive coal outcrops occurring in the sedimentary rocks of the Lignitifero Formation (Fig. 1). During this period, coal became the primary national resource and was mostly exploited in open-pits (Fig. 2).

Location map (left); Geological sketch map (right); The geological cross-section (below)

Location of exploitation panels, W1, W2, W3, W4. On right the maps of underground tunnels and excavation areas

Since 1976, Carbosulcis S.p.A., which is the current owner of the mining licenses in the Sulcis district, started mining activity south of Seruci, in the Monte Sinni area, nearby the Nuraxi Figus village (Fadda et al. 1994). The Monte Sinni coal panels are located at a depth of about 400 m b.g.l., extend for an area of 55 km² and are mined through 30 km of tunnels (15 km of which correspond to permanent infrastructures). Mining activity was conducted in correspondence with some exploitation panels, called W1, W2, W3, etc., with average dimensions of 300 m × 600 m (width × length) and a height of about 3 m (Fig. 3). In the initial period of activity, mining was developed through room and pillars method, whereas from the 1980s the coal extraction was carried out with a different method, called “longwall cutting in retreat”. The latter is based on the continuous extraction of the coal panel through a shearer that digs along the coal face, while a transporter belt carries away the grained coal. When the shearer moves on, the roof of the gallery at the back of the active face collapses. The Carbosulcis S.p.A. used the longwall method to extract the W3 and W4 panels, in the periods between 2008 and 2010 and from 2011 to 2012, respectively. In 2018, following the progress of environmental regulations on the cessation of the use of coal, the Company started the closure plan that progressively will bring to reclaim all the areas by 2027.

Flow chart of the methodological approach

Geology of the Sulcis district is characterized by Cenozoic sedimentary and volcanic rocks, unconformably overlaying a Paleozoic basement. The Cenozoic sedimentary rocks consist of limestones, sandstones, conglomerates, marls, and silty clays, and are formally subdivided into four stratigraphic formations. At the basis of the Cenozoic succession, the oldest rocks are represented by the limestones of the Macro-foraminifera Formation, which are covered by the sandstones, marls and limestones of the Miliolitico Formation (20–70 m thick; Lower Eocene age). The following stratigraphic interval is represented by clays, marly limestones, bituminous limestones, marls, and conglomerates of the Lignitifero Formation (70–150 m thick; Lower-Middle Eocene age). which hosts the coal horizons object of mining activity. Coal seams commonly are 10 cm thick, rarely reaching 30–50 cm, and are interbedded with clays (Assorgia et al. 1992a and b). The Lignitifero Formation is covered by about 300 m of sandstones, conglomerates, and marls of the Cixerri Formation (Eocene–Oligocene; Pasci et al. 2010). The whole Cenozoic sedimentary succession is in turn covered by volcanic rocks, represented by twelve andesitic, dacitic to rhyolitic ignimbrite plateaus belonging to several stages of explosive volcanic activity of Oligo-Miocene age (Morra et al. 1994). The whole district has been dislocated by normal faults, which lower the sedimentary succession to various hundred meters of depths within the basin (Barca et al. 2000). These geological features allow the existence in the area of two types of aquifers: one occurring in the volcanic rocks and a second one in the Miliolitico Formation, which have been both dewatered by Carbosulcis S.p.A. before the start of the underground mining activity.

From a geotechnical point of view, the Rock Mass Rating (RMR) (Bieniawski 1989) was used to numerically define the characteristics of the rocks, using "RMR value" for defining five quality classes (from very poor to very good). In particular, 89% of the volcanic rocks correspond to the fair class (III), whereas the sedimentary Formations are characterized by 57% to the fair class (III), 28% to the poor class (IV), and 10% to the good class (II) (Fadda et al. 1994). According to Fadda (1994) the Cixerri Formation has a heterogeneous structure, where conglomerate and marls layers belong to the lowest classes, while sandstones are characterized by the best performance.

For the present study, as summarized in the flow chart of Fig. 3, satellite data and mining information were used to create a near-real time monitoring system. Satellite data were collected in ascending and descending geometries covering the time span 2011–2020, to study subsidence by the vertical deformation maps and the displacement time series. Two datasets composed of high-resolution images have been used:

1. (1)

   Satellite data 1 (SD1): PSP-IFSAR algorithm (Costantini et al. 2008), time 2011–2014, obtained in the framework of Not-Ordinary Plan of Environmental Remote Sensing (Piano Straordinario di Telerilevamento Ambientale—in Italian), funded by the Italian Ministry of Environment (Costantini et al. 2017; Di Martire et al. 2017).

2. (2)

   Satellite data 2 (SD2): SUBSIDENCE software (Mora et al. 2003; Iglesias et al. 2015), covering the period 2013–2020.

SD1 consist of 41 images acquired in ascending geometry and 58 images acquired in descending geometry in the time interval May 2011–March 2014. SD2 contains 102 and 116 images in ascending and descending orbit, respectively related to the period from October 2013 to July 2020. First data set derived from the interferometric processing within the framework of the third stage of the PST-A project (Not-Ordinary Plan of Environmental Remote Sensing) funded by the Ministry for the Environment and Protection of the Territory and the Sea (Italian National Geoportal, www.pcn.minambiente.it). The second dataset was obtained by the processing of SD2 using SUBSIDENCE software, which implements the Coherent Pixels Technique (CPT) algorithm, developed at the Remote Sensing Laboratory (RSLab) of the Universitat Politecnica de Catalunya of Barcelona. SUBSIDENCE uses the CPT-Temporal Phase Coherence (CPT-TPC) approach to extract from a stack of differential interferograms the deformation evolution over wide areas during large periods (Mora et al. 2003; Iglesias et al. 2015). The processing is structured in three main phases:

(1) Interferogram generation: the generation of the best interferogram set among all the available images of the zone under study;

(2) The Stable Coherence Scatterers (SCS) selection: the points in the detected area characterized by signal stability and higher than the threshold in a specific percentage of interferograms. To get enough SCSs a coherence limit of 0.6 was considered, assuming an error in mean displacement rate lower than 1.5 mm;

(3) Linear velocity of deformation (Iglesias et al. 2015): linear deformation time series (TS) were calculated starting from phase analysis. The Delaunay triangulation was used to check an offset among the different interferograms. Then, it was applied a filtering process to assess the deformation evolution of selected pixels (SCS) (Blanco et al. 2008); the geocoding of the results was carried out in WGS84-UTM.

This process allowed to elaborate PS (Permanent Scatters) maps for each acquisition geometry in terms of mean displacement rate and time series of deformation, along the Line of Sight (LoS) of the satellite (Fig. 4). For a better understanding of the results, it is necessary to clarify the acquisition mode. Satellites move two semi-polar orbits defined as ascending and descending, moving from South to North and North to South, respectively, sending electromagnetic signals along the Line of Sight (LoS) to identify targets on the surface. By convention, displacements approaching from the sensor will take on positive values (blue), while displacements away the sensor will take on negative values (red—Fig. 4).

Geometry acquisition: a Ascending; b Descending

After this, by using images acquired in ascending and descending geometries the vertical displacement component was calculated (Cascini et al. 2010; Di Martire et al. 2013). Subsequently, to improve the analysis of the Nuraxi Figus subsidence, satellite datasets were integrated with the following information, provided by Carbosulcis S.p.A.: stratigraphic logs of 9 drillholes, geotechnical parameters of country rocks, mining reports regarding the exploitation of panels W3 and W4, and maps of underground tunnels and excavation areas. Spatial migration in time of the excavation front within the mining panels W3 and W4 was compared with the evolution of vertical ground displacements detected with the satellite data during the mining activity.

The 3D geological model of the subsurface was elaborated by using lithostratigraphic data from 7 boreholes through the software Rockworks® (Fig. 5).

a Location of stratigraphic logs used to elaborate b 3D geological model of subsurface

The model allowed us to determine the thickness of the overburden strata in the various parts of the study area. A stratigraphic section perpendicular to the W3 and W4 exploitation panels was extracted from the 3D geological model to produce a numerical subsoil model through the Plaxis® 2D software (Brinkgreve et al. 2008). The finite element numerical model allowed us to model the excavation process in two-dimensions in order to quantify the subsidence induced at surface level. The characteristic values of geotechnical rock parameters used for modelling (Table 1) were obtained by Carbosulcis S.p.A.

As reported in Table 1, due to similar geotechnical parameters, the twelve ignimbrite plateaus overlying the sedimentary units were grouped in four geotechnical layers. Cross sections (width = 200 m, height = 3 m) of panels W3 and W4 were considered in the 2D model at the depth of 400 m. Each rock layer has been considered in the numerical model as linear elastic, with Young modulus E and Poisson’s ratio n as reported in Table 1. Evaluation of the overall subsidence (total vertical displacements of the surface) has been performed at the end of the following calculation steps: (1) initial in situ stresses (geostatic stresses), (2) excavation of W4 panel, and (3) excavation of W3 panel.

5.1 Processing phase

Starting from two satellite datasets, four PS-displacement rate maps for both acquisition geometries (ascending, descending) have been created. The PSs obtained were imported into GIS platform and subsequently made visible according to a colour scale: negative values conventionally indicate a movement of the target away from the satellite (LoS), while positive values indicate movement towards the sensor; stable areas are shown by using the green color (see Fig. 4). The maps represent the average velocity recorded during the period covered by the acquisitions (Fig. 6). The results of both ascending and descending geometries show that the predominant displacement rates are located in correspondence with the panels.

Mean displacement rate maps in ascending and descending orbit

In the central part of the monitored area, in the first period considered (2011–2014), PS-deformation velocities exceed the detection limits (Colesanti et al. 2006), and no satellite-based measurements are available. On the contrary, in the second period analyzed (2013–2020), the same zone is covered by PS-measurements.

This fact is likely related to the decrease in displacement velocities. Indeed, the cumulated satellite-based LoS displacements vary between − 130 and + 28 mm (until January 2014), − 293 and + 28.4 mm (until March 2014) and between −6.9 and + 1.6 mm (until March 2020), − 8.72 and + 4.33 mm (until March 2020) in ascending and descending geometries, respectively. In the 3D subsurface geological model of the investigated area (Fig. 5), it is possible to see the position of the Lignitifero Formation, which hosts the mined horizon at 350–500 m b.g.l. The thickest rock layers overlying this interval correspond to the sandstones and conglomerates of the Cixerri Formation. The volcano-pyroclastic rocks and ignimbrites (named Volcanic 1 to 4) have a thickness that varies from 100 to 298 m.

5.2 Post-processing phase

The numerical model produced a simulation of the expected subsidence induced at the surface by the collapse of voids created by the excavation of the two coal panels (W3 and W4). The model was produced along an N-S section crossing both the W3 and W4 panels (Fig. 5). To reproduce the chronology of the events and analyze if the differential exploitation of the two panels influenced the total measured subsidence, two calculation steps were performed. The first was related to the sole exploitation of the panel W4 (occurred from 2008 to 2010), whereas the second was related to the excavation of the W3 panel (which occurred from 2010 to 2012). Figure 7 shows that vertical displacements reach a maximum value of 32 cm. The lateral extent of the subsidence covers a horizontal length of about 1100 m, which is quite wider than the mined panel size. At the surface, the expected vertical displacement reaches a maximum value of 18 cm.

Vertical displacements obtained by means of the geotechnical numerical model. On top subsidence profile related to the section AA’

The availability of both ascending and descending datasets allowed us to reproduce the kinematics of the deformation through the displacement vector decomposition. The results, as shown in Figs. 8 and 9, are maps of vertical (VC) and horizontal (HC) components, with the total amount of displacement calculated in the time span 2011–2020. Considering an x, y, z Cartesian coordinate system, HC maps coincide with the horizontal W-E and the VC maps with the vertical components. VC and HC were calculated considering three phases: pre-mining, syn-mining and post-mining activity. The availability of satellite dataset allowed to obtain one map during the pre-mining activity related to September 2011, where no displacement has been detected. During the W3 panel extraction carried out in the period 2011–2012 the displacement vector components were analyzed every two months. To understand the temporal evolution of displacements the progress of mining exploitation (i.e. the position of the exploitation front) has been integrated into the vector maps.

Vertical deformation maps pre-mining, syn-mining, and post-mining activity. Black dashed lines in the figure represent the mining panels

Horizontal deformation maps pre-mining, syn-mining, and post-mining activity. Black dashed lines in the figure represent the mining panels

The post-mining activity was investigated every year from January 2013 onward. It is important to note that during the excavation of the W3 panel, between November 2011 and October 2012, only slight deformations can be identified. In detail, displacements started in November 2011 and developed until October 2012 with maximum cumulated vertical displacements of ca. 11 cm and horizontal ones between −14 cm and 4 cm. On the contrary, the largest deformations were detected in correspondence of post-mining operations, during the period from January 2013 to January 2014. Subsequently, subsidence continued slowly until January 2020 with maximum cumulated vertical displacements of ca. 26 cm and horizontal ones between − 24 cm and 7 cm.

The aim of this work was to create a near-real time system for monitoring mining activity starting from PS-maps and using a subsurface3D geological model. The expected subsidence value was defined in the geotechnical numerical model of the rocks overlaying the two mining panels W4 and W3. Subsequently, the temporal vertical and horizontal component maps were created for monitoring the mining area during the exploitation and post-exploitation activity. Figure 10 shows the temporal evolution of vertical displacements along the N-S oriented profile cross-cutting the mined panels, compared with the expected subsidence along with the same profile, as it was generated by the numerical model. The deformation started in November 2011, in correspondence with the W4 panel, and developed and extended until 2020, with maximum cumulated vertical displacements of ca. 26 cm located in the centre of two panels. It is important to note that most of the subsidence was temporally delayed respect to the coal exploitation.

Cumulated vertical deformation profiles during the pre-mining, syn-mining, and post-mining activity. Yellow dashed lines in the figure represent the expected subsidence

In fact, during the mining activity, the expected subsidence was not exceeded, whereas the maximum deformation occurred during the post-mining operations. The delayed occurrence of the maximum deformation depends on several factors: peculiar features of the overburden geology, thickness of the panels, depth of the excavation area, and characteristics of the mining method. Several coal mining districts located in different countries, characterized by variable geological settings, show temporal delays in the occurrence of maximum surface deformations. This particular subsidence phenomenon, which is called residual subsidence (Alheib et al. 2005), continues after the extraction and in some cases can occur some months or years after the end of the underground mining activity (Huang et al. 2020; Cui et al. 2020; Modeste et al. 2021). Sometimes the residual surface deformation could be also characterized by uplift phenomena (Vervoort 2020). To understand the surface displacement, in Table 2 maximum horizontal (E-W direction) and vertical components are compared with mining progress. During the extraction in 2011–2012, the major displacement reported showed a horizontal component (Fig. 9). This was likely due to the position of the new mined zone. Subsequently, it is possible to note that the vertical displacement increases respect to the horizontal one (Fig. 8). After the conclusion of mining activity, the deformation is concentrated between the two panels in the area characterized by the greatest geotechnical weakness (Table 2). This confirms that the observed displacements are directly related to the excavation process. The W3 panel was mined from 2011 to 2012 and the highest deformation increment was detected in January 2013, with differential vertical displacements of ca. 10 cm. In later times, the increment is of a few centimeters per year, with values decreasing with time until now.

According to several authors (Cui et al. 2000, 2020; Alheib et al. 2005; Tajduś et al. 2021), subsidence due to longwall mining is structured into three intervals. The first called initial subsidence, is the period when the surface movement lightly starts with low velocity.

The second phase (principal) is considered when displacements increase, generally up to 80% to 90% of final subsidence. The final phase, named “delayed” or “residual” period, is associated with surface subsidence that continues after the end of mining activity, and can be characterized by 10% to 15% of final subsidence. This “delayed” period can start around 12–18 months, or also around 3–4 months after the end of the underground excavation, depending on the high or low geotechnical quality of the overburdened strata, respectively. In the Nuraxi-Figus study case, the second subsidence phase ends one year after the end extraction (Fig. 11), when the vertical displacement achieves about 84% of the final subsidence. The delayed subsidence starts in January 2015 and develops slowly with a rate of about 20 mm/a (7% of final subsidence).

Profile of maximum vertical cumulated displacement compared to three subsidence phases: initial (green), principal (red), residual (yellow)

In the initial stages of underground mining, the caving zone is characterized by elastic deformations that recede if the overburdened forces are deleted. If extraction continues after a certain threshold, deformations become permanent. In this case, the pressure on the exploited areas increases, the sides move inward, the floor is subjected to uplift and the roof slides down (National Coal Board book 1965; Shadbolt et al.1978). According to Whittaker (1989), deformation occurring above a mining goaf depends on the following factors: depth of the cover, properties of overburden strata, seam thickness, geometry of the extraction panel, surficial topography, and extraction techniques. Specifically, in the underground mining four zones can be identified in the rock pile occurring above the extraction area. Starting from the bottom to the top there are: I—the zone immediately above the extraction area (caving zone), II—the fractured zone where the major cracks are present; III—the deformation zone, and IV—the surficial zone where subsidence could occur (Peng 1992). According to Mills (1998), the maximum subsidence (Sm) expected in an underground mining area is related to the width of the extraction panel (W), and height of cover strata (H). By using these parameters, it is possible to identify three types of excavation areas (Whittaker et al. 1989; Mills 2009):

1. (1)

   Supercritical (W/H > 1.6): Sm could be between 55% and 65% of the mined seam thickness;

2. (2)

   Critical W/H (from 0.6 to 1.6): Sm could be about 10% of the mined seam thickness, depending on changes in panel geometries, the overburden depth, and the composition and geotechnical properties of the strata;

3. (3)

   Subcritical (W/H < 0.6): Sm is negligible.

Applying the above subsidence model to the study area, it is possible to say that the caving zone should be completely included in the Lignitifero formation, whereas the fractured zone should affect the lower part of the Cixerri formation. The residual overlaying part of the Cixerri formation and the volcanic rocks should belong to the bending zone. In the geotechnical numerical model, it is possible to note that the maximum displacement is related to the Cixerri Formation. This is probably due to pre-existing structural discontinuities (fractures and faults), occurring in this Formation (Fadda et al. 1994), that were not considered in the geotechnical model. The surface zone is instead characterized by slight surface movements.

In this study, an integrated approach for near-real time mining monitoring was presented. The method allowed us to detect a subsidence phenomenon that occurred during the mining activity and continued also for a few years after the completion of mining work. The study confirms that ground surface deformations can occur also in areas subjected to very deep underground mining and that the DinSAR techniques produce reliable results in terms of monitoring in a certain time span. The geotechnical numerical method can be used as subsidence forecasting model allowing to obtain the maximum expected amount of ground displacement. In underground mining areas, it is important to use monitoring techniques that can provide cost-effective regional perspectives for preventing possible negative environmental effects. The near real-time monitoring system applied in this study proved to be very useful for detecting subsidence during the mining activity and in the post-mining period. It is also important to emphasise that, in the absence of on-site monitoring, the method applied allowed us to reconstruct the phenomenon that occurred in order to undertake the necessary activities for the rehabilitation of the area. Therefore, it can be used for planning environmental remediation plans. However, to support the previous analysis and better investigate the subsidence-related cause-effect mechanism, as effectively as possible, it might be necessary to also acquire additional external data (i.e. leveling surveys, extensometers).