Monitoring Land Surface Deformation Associated with Gold Artisanal Mining in the Zaruma City (Ecuador)

Abstract

Underground mining can produce subsidence phenomena, especially if orebodies are surficial or occur in soft rocks. In some countries, illegal mining is a big problem for environmental, social and economic reasons. However, when unauthorized excavation is conducted underground, it is even more dangerous because it can produce unexpected surficial collapses in areas not adequately monitored. For this reason, it is important to find quick and economic techniques able to give information about the spatial and temporal development of uncontrolled underground activities in order to improve the risk management. In this work, the differential interferometric synthetic aperture radar (DInSAR) technique, implemented in the SUBSOFT software, has been used to study terrain deformation related to illegal artisanal mining in Ecuador. The study area is located in Zaruma (southeast of El Oro province), a remarkable site for Ecuadorian cultural heritage where, at the beginning of the 2017, a local school collapsed, due to sinkhole phenomena that occurred around the historical center. The school, named “Inmaculada Fe y Alegria”, was located in an area where mining activity was forbidden. For this study, the surface deformations that occurred in the Zaruma area from 2015 to 2019 were detected by using the Sentinel-1 data derived from the Europe Space Agency of the Copernicus Program. Deformations of the order of five centimeters were revealed both in correspondence of known exploitation tunnels, but also in areas where the presence of tunnels had not been verified. In conclusion, this study allowed to detect land surface movements related to underground mining activity, confirming that the DInSAR technique can be applied for monitoring mining-related subsidence.

1. Introduction

In several countries, mining activity is vital for social and economic development. Sometimes, the exploitation of deposits is related to resources of high economic value, increasing the interest both of mining companies and artisanal miners. Artisanal and small-scale mining (AMS) is recognized as a considerable source of revenue for millions of people in about 80 countries worldwide [1]. ASM takes place in several regions of the world, mostly in the global South—sub-Saharan Africa, Asia, Oceania, Central and South America. In 2017, about 40 million people worldwide were working in the artisanal mining sector and the industry had a significant impact on global mining [2]. These small-scale operations, which are typically undertaken outside the law and without compliance with industry standards, can promise great rewards. However, violations of rules are quite common, and these activities can damage and pollute the territory environment. In certain environments, mining activity can produce subsidence phenomena [3,4] and can be even more dangerous if it is carried out illegally because it can produce unexpected sinkholes [5,6,7] in areas not adequately monitored. Monitoring and parameterization of these dynamic processes are particularly relevant. Over the last years, several monitoring methods have been used, including remote sensing techniques. Remote sensing has proved to be highly suitable for these and for many other necessities, since it captures data covering vast study areas. Among the different types of remote sensing techniques, one that very well meets dynamic processes requirements is the differential interferometric synthetic aperture radar (DInSAR). In recent years, DInSAR techniques have been effectively used to measure surficial deformations related to subsidence phenomena [8] and to monitor underground mining activity [9,10,11,12,13]. The hills of Zaruma and Portovelo (El Oro province, Ecuador) have been mined for gold and silver for centuries. The Incas were already extracting gold and silver in the area with hydraulic mining of the oxidized parts of veins. In 1549, Mercadillo, one of Pizzarro’s force, following upstream the Amarillo River encountered the Inca mine and founded the town of Zaruma [14]. Since then, legal and illegal mining has been the reality of the area, bringing many problems that range from poor land and city planification, to criminal and illegal activities. In the last 10 to 15 years, a big concern in the city and its surroundings is represented by land subsidence related to the combination of geology, morphology and mining activity [15]. In fact, one of the principal triggers of this phenomenon is considered the presence of old/new legal/illegal galleries, shafts and other underground excavations. In this work, we have studied surficial deformations that occurred in the Zaruma historical center, between 2015 and 2019, in order to compare such deformations, on one hand, to known tunnel mining and, on the other, to detect possible locations of unknown tunnels. To pursue this aim, we used SENTINEL-1 image processing derived from the Europe Space Agency of the Copernicus Program.

The work is structured as follows: first the mining activity at the Zaruma site is described, then a geological setting of the area is performed. A brief description of the data and the interferometric technique are carried out. Finally, the results are reported, and conclusions are drawn on the application made and possible future developments.

2. Zaruma Geological Setting and Mining Activity

The city of Zaruma or “Villa del Cerro del Oro de San Antonio de Zaruma” is located in southwestern Ecuador. The geology of the study area is characterized by pre-Mesozoic metamorphic rocks, which comprise the Amotape-Tahuin (Ta) massif in the south, and the Chaucha metamorphic and volcanic rocks in the north [16] overlaid by younger volcanic rocks. In the Zaruma area, the basement, called the El Oro metamorphic complex, is covered by calc-alkaline volcanic materials of the Saraguro formation (FmSa) (Figure 1a,b). The FmSa is widespread all over southern Ecuador, characterized by variable lithology, comprising andesitic and andesitic basalt lavas, highly weathered andesitic tuffs and dacitic tuffs, and ignimbrites [17,18].

The pyroclastic rocks range from fine-grained tobaceous to agglomerated (tuffs and thick agglomerates with blocks of lava) and alternating porphyritic andesite lavas [18,19,20]. Another lithology occurring within the study area is represented by the Zaruma Urcu Rhryolite (RZU), mostly consisting of dyke stocks and other forms of rhyolitic intrusions [21]. Small occurrences of metamorphosed gray sandstones, interspersed with limolites and dark gray shales, likely corresponding to the Paleozoic Ta, can be locally observed [17,19,20]. According to the results of geological and geotechnical studies, the predominant volcanic rocks in the study area are characterized by a high weathering degree, very intense surficial reworking ranging from poor to very poor quality. These characteristics notably increase terrain instability. The epithermal gold orebodies mined at Zaruma are genetically associated with calc-alkaline Eocene to Late Miocene igneous complexes occurring in central and southwestern Ecuador [14,21]. Mineralization within the area is considered to be an intermediate to low-sulphidation epithermal to a mesothermal gold-silver-lead-zinc-copper system [22,23]. In the Portovelo-Zaruma area, the gold-bearing quartz vein system shows a north–south trend, characterized by sub-parallel structures exclusively located within the Cretaceous semi-high altered andesitic rocks. Three main types of gold-bearing veins are present in the study area (Figure 1c): the quartz-pyrite zone: quartz veins with disseminated pyrite, minor chlorite strikes, bands and patches; the sulphide zone: quartz veins with high pyrite, chalcopyrite, galena and sphalerite in bands, patches and coarse disseminations; and the calcite zone: carbonate veins with coarse calcite and calcite-quartz, galena, sphalerite and chlorite in occasional nodules [14,17,18,19]. These metallogenetic features make the Zaruma area enhanced with several gold mineralization. The city is well known because of gold mining as most of the population works or lives nearby and is in contact with mining zones [24]. The mining activity of Zaruma started during the 16th century because of the presence of gold veins and big quantities of gold in the rivers, and the city took the name of “Asiento de Minas de Zaruma”. Recently, the mining industry in Zaruma has gone through different phases. In the early 20th century, the mining activity was first controlled by the company SADCO and then by CIMA corporation, but also by individuals carrying out artisanal activity, following the main gold vein called “Sweet Water” [25]. Later on, a lot of people working in the gold exploration and exploitation in Zaruma independently built their houses, causing a disordered and accelerated growth of the urban area [24]. Contextually to the licensed mining activity, widespread illegal gold exploitation all around the city area has been developed and many tunnels and galleries were built there. The uncontrolled mining increased the environmental and economic dangers in Zaruma. One of the most relevant infrastructure problems that occurred in Zaruma is the collapse of the school, named “Inmaculada Fe y Alegria”, in 2017.

The school collapsed into the ground after a big sinkhole episode (Figure 2a), and over the years other problems have occurred in the surrounding areas (Figure 2b,c). Most of the mentioned problems were related to illegal artisanal mining, carried out with low and poor technical standards and minimum control, along the old SADCO and CIMA galleries, which were built without any regulation.

In order to manage problems of illegal mining that affect the cities of the El Oro province, the Ministry of Energy and Non-Renewable Natural Resources and the Risk Management Services (SGR) organized a plan called “Zaruma - Portovelo 2022” [29]. The project has several actions aimed at studying the deformation due to outlaw activity and restore the urban area. Among the first actions implemented by the SGR, geoelectrical surveys were carried out in the Zaruma historical center, which allowed a preliminary identification of some tunnels. Subsequently, other studies, carried out by the University of Cuenca Politecnica Salesiana, allowed to obtain an updated map of the tunnels [30] (Figure 3). It is worth to point out that these studies were certainly not exhaustive, as they were only limited to the historical part of the city.

3. Materials and Methods

It is now almost three decades (early 1990s) that the differential interferometry SAR technique (DInSAR) [31] used for monitoring surface deformations has caught the attention of the scientific community. DInSAR, in fact, is a very effective technique for the accurate measurement of slow ground movements, structures and infrastructures due to subsidence, landslides, earthquakes and volcanic phenomena [32,33,34,35]. By means of DInSAR, it is possible to measure ground displacements with sub-centimetric accuracy, starting from data acquired from satellites orbiting the earth at an average height of 600 km. The DInSAR approach is based on the analysis of phase difference in interferometric stacks of radar images: this technique operates at a full spatial resolution and identifies reliable scatterers (permanent scatterers—[36]) by measuring their multitemporal coherence related to the phase stability. In particular, DInSAR techniques allow to analyze long data series producing mean displacement rate maps and time series of deformations along the direction between the SAR sensor and the target (line of sights—LoS). In this work, TOPS IW Sentinel-1A and B ascending and descending mode images, obtained via the Sentinel Scientific Data Hub [37], acquired in the time span from June 2015 to June 2019 (

Table 1. Synthetic aperture radar (SAR) data stacks analyzed in this study.

Satellite	Orbit	Period	Nr Scenes
Sentinel1 A/B	Ascending	17/01/2016–12/062019	89
Sentinel1 A/B	Descending	05/062015–23/092018	51

), were processed by the SUBSOFT software, which uses the coherent pixel technique-temporal sublook coherence (CPT-TSC) approach [38,39], and developed at the Remote Sensing Laboratory (RSLab) of the Universitat Politecnica de Catalunya of Barcelona.

Specifically, the interferometric chain implemented in SUBSOFT is divided into the following main steps: the first is the generation and selection of interferograms, which consists of selecting a quality set of interferograms, considering the temporal and spatial thresholds, from all the possible combinations among available images. The second step is the pixel selection, named stable coherence scatterers (SCS—[39]), characterized by a reliable phase. In order to select SCS, a TSC estimator [39] has been used.

The third one is the evaluation of the linear term of deformation to define the linear velocity and location of the targets; then the evaluation of the non-linear deformation component and atmospheric artifacts to assess the deformation evolution of selected pixels; and finally, the geocoding of the results in WGS84 and WGS84-UTM. Further details can be found in [38,39].

After the processing step, the interferometric results were post-processed by the application of the kernel density estimation (KDE) algorithm [40], which allowed to identify in a rapid way the unstable areas (UAs) affected by meaningful deformations [41,42,43,44,45].

Specifically, the “Heatmap” plugin in the Qgis software has been used. This creates a density map (heatmap) of an input point, using the following formula:

$$\mathrm{KDE}=\frac{1}{{\left(Bandwidth\right)}^2}{\displaystyle \sum }_{i=1}^n\left[\frac{3}{\pi }\left(w_i\right){\left(1-{\left(\frac{L_i}{Bandwidth}\right)}^2\right)}^2\right]$$

(1)

$$for L_i<Bandwidth$$

(2)

where

• Bandwidth is the radius;
• w_i is the weight;
• L_i is the distance between point i and the (x, y) location.

The density is a function of the number of points, or the sum of the weight parameter. In fact, the use of this parameter allows to increase the influence that some points have on the resulting density map. In this study, the mean displacement rate is considered as a weight parameter.

Heatmaps allow quick identification of “hotspots” of points. As for the kernel map, the first feature to set is the choice of bandwidth or the search radius. The bandwidth controls the smoothing of the results, larger values result in greater density, but smaller values may show finer details. Therefore, it is important to estimate the bandwidth starting from the presence data available in order to use the best radius value.

In order to calculate the best radius that must be used in the kernel density, the formula of [40] has been adopted:

$$Bandwidth=0.9\times \min (\mathrm{SD},\sqrt{\frac{1}{\ln \left(2\right)}}\times D\mathrm{m})\times n^{-0.2}$$

(3)

$$*\mathrm{SD}=\sqrt{\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i}}{\left({\mathrm{x}}_{\mathrm{i}}-{\mathrm{X}}_{\mathrm{w}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i}}}+\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i}}{\left({\mathrm{y}}_{\mathrm{i}}-{\mathrm{Y}}_{\mathrm{w}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i}}}}$$

(4)

where

• Min is the minimum value between the two in the interlude;
• SD is the standard distance;
• D_m is the (weighted) median distance from (weighted) the mean center;
• n is the number of points;
• x_i, y_i are the coordinates;
• w_i is the weight.

The radius was calculated using Qgis tools: in detail, the centroid was considered through the Mean Coordinate tool and then the Raster Calculator was used to estimate the bandwidth for each satellite geometry. The values obtained were inserted in the Heatmap tool and the mean displacement rate as the weight field has been used (w_i) to increase the quality of the results.

Finally, thanks to the availability of images acquired in the two ascending and descending geometries, it has been possible to obtain the vertical (Z) and E–W components of the displacements. This operation is possible (system resolution of 2 equations in 2 unknowns (Equation (5))) under the hypothesis that the motion component in the north–south direction is negligible [46,47]. This is an intrinsic limit of the acquisition system, since the north–south motion component is strongly underestimated since the SAR sensor’s direction of view is almost orthogonal to the north–south direction (the satellite travels along quasi-polar orbits with antennae oriented orthogonally to the flight direction) and is therefore not very sensitive to variations in this direction.

$$\{\begin{array}{c}{D}_a=D_x{s}_{xasc}\times D_y{s}_{yasc}\times D_z{s}_{zasc}\\ D_d=D_x{s}_{xdesc}\times D_y{s}_{ydesc }\times D_z{s}_{zdesc}\end{array}$$

(5)

where D_a and D_d are, respectively, the displacement values in ascending and descending geometry, D_h and D_v are the displacement vector components along the horizontal (E–W) and vertical directions and S_xasc, S_yasc, S_zasc, S_xdesc, S_ydesc, and S_zdesc are the incidence angles in the two geometries.

4. Results

The datasets available from this study consist of 89 images acquired in ascending geometry with a time-revisiting variable between 12 and 36 days in the time span January 2016–June 2019, and 51 images acquired in descending geometry with a time-revisiting variable between 12 and 96 days in the time interval June 2015–September 2018. From all possible pairs of interferograms, only those characterized by spatial and temporal baseline thresholds of 50 m and 150 days, respectively, have been selected. In particular, 793 and 271 interferograms in ascending and descending orbits, respectively, have been identified (Figure 4).

The following step consisted of SCS selection. CPT-TSC allowed to select points in the detected area characterized by a phase quality higher than a threshold value set by the operator according to the error in the displacement evaluation considered acceptable (in this case less than 1.5 mm), which in turn is a function of the expected mean displacement rate. In this case, a TSC value equal to 0.65 has been set in both geometries (Figure 5) in order to obtain an acceptable displacement error, lower than 1.5 mm, and to select an adequate number of points.

Using the CPT-TSC approach, two mean displacement rate maps have been obtained, each for both acquisition geometries. In the El Oro province territory about twenty and twenty-one thousand targets (ascending and descending) located in the Portovelo, Piñas and Zaruma districts have been detected. Specifically, in the Zaruma town, about 4000 and 5200 targets have been detected, in ascending and descending orbits, respectively. The results allowed to evaluate the surface ground displacements in terms of mean displacement rate and time series of deformation along the LoS of the satellite. Figure 6 shows the mean displacement rate in terms of cm\year recorded in the time span 2015–2019. The maps have been represented using a color scale from red to blue, where the negative values conventionally indicate a movement of the target towards the satellite, while the positive values indicate movement far from the sensor, and the green color identifies stable areas.

In order to identify the displacement rate threshold for which areas should be considered stable, a coefficient of variation, given by the ratio between the standard deviation and the velocity average module, for all selected SCS, has been evaluated. The threshold value obtained is about 0.3 cm/year, corresponding to a value for which standard deviation is higher than the mean displacement rate value [46].

To detect the presence of UAs in the urban area of Zaruma, the KDE algorithm was applied to both analyzed datasets (Figure 7). Such a figure shows the hotmaps in the area of the Zaruma historical center displaying the density of targets weighted on mean displacement rate values during the covering period for each satellite orbit. Both blue and red hotmaps highlight where surface deformations are present. As already described above, the color blue indicates movements far from the satellite, while red movements are towards the satellite. Moreover, the intensity of the color gives information on the magnitude of the deformations, while the radius indicates the potential extension of the phenomenon.

In detail, the detected UAs are distributed as reported in

Table 2. List of unstable areas (UAs) and probable cause of the observed SAR displacement rates.

Unstable Area	Probable Cause
UA1	Reported tunnels
UA2	Reported tunnels
UA3	Reported tunnels
UA4	Unknown tunnels/Landslide
UA5	Unknown tunnels/Landslide
UA6	Unknown tunnels
UA7	Reported and unknown tunnels
UA8	Unknown tunnels
UA9	Unknown tunnels
UA10	Unknown tunnels

: (a) south of the Inmaculada school (UA1), where displacement rates of about 3 cm/year have been recorded and where tunnels are known; (b) around the school (UA2), which collapsed in 2017, where significant deformation rates can still be observed (about 1 cm/year). It is worth to point out that probably higher displacements affected the school during the 2017 event. Such displacements were not detected, related to their extreme rapidity in a short time span, not allowing the interferometric technique to detect a coherent signal because of its inability to analyze such rapid deformations [46]. Meanwhile, it was possible to verify that the subsidence phenomenon affecting the area was already in place before the collapse; (c) on the eastern hill of the inhabited center (Mirador La Colina—UA3), with mean displacement rates of the order of 1 cm/year, where tunnels have been mapped; (d) UA4 has been identified to the south, in correspondence of the Humberto Molina Hospital with rates of the order of 1 cm/year. In this case, due to the absence of known tunnels and the morphology of the site, these displacement rates could be caused by landslides (rotational slide) which affect the slope or subsidence induced by the presence of lawless activities; (e) as for UA4 and the UA5, located on the west side of the Marcelo Zambrano Stadium, they can be associated with a landslide phenomenon (rotational slide), due to the presence of a steepness slope or lawless activities; (f) UA6, located in correspondence of the Marcelo Zambrano Stadium, is characterized by null slopes, and by the absence of reported tunnels. Therefore, these peculiarities lead to the interpretation of recorded deformation rates due to the presence of illegal activities; and finally, (g) UAs 7 and 10, in the northern sector of the town and UAs 8 and 9, in the westernmost sector of Zaruma (near the IESS Hospital), showed mean displacement rates of the order of 1–2 cm/year, even though there are no tunnels in the official inventory. It was not possible to validate the results by comparing them with in situ measurements, due to the absence of the latter, but only by consulting several newspaper articles and scientific papers [15,26,27,28,30].

Finally, a comparison between UAs and reported mining tunnels have been carried out, highlighting that out among the 10 UAs identified, 40% of these have at least one tunnel mapped in the official inventory, while for the remaining 60% there are no reported tunnels, so the presence of lawless activities could be assumed.

5. Discussion

The results obtained in this work have, on the one hand, confirmed the phenomena of instability, induced by the presence of underground activities, legal and illegal, which affect a large part of the urban center of Zaruma, and on the other, highlighted the need to implement a deformation monitoring system in order to mitigate the risks present. Actually, the aim of this paper was to evaluate the capability of Sentinel-1 data to detect and characterize the possible deformations due to mining activity. The processing of interferometric data by means of the DInSAR technique and the subsequent interpretation (kernel density) allowed to identify unstable areas (UAs) affected by subsidence and other gravitational phenomena. The use of the interferometric technique to monitor mining areas affected by subsidence phenomena which induce damage to structures has been carried out in several works. Many papers focus on areas where mining activities are known and often where in situ monitoring systems are implemented [12,13]. The most complex situation is that of unauthorized mining areas. In this case, the main purpose is the identification of areas potentially subject to subsidence and then the deformation rate. An example is a recent work [15], implemented in Zaruma, where the authors, by processing a limited number of images only in ascending orbit, have identified potential unstable areas, for an extension of about 40% of the urban center, but without evaluating their temporal evolution. Unlike what [15] carried out, in this work, the images in both acquisition geometries (ascending and descending) were processed for a longer time interval (2015–2019). Moreover, a post-processing phase, using the kernel density algorithm, has been implemented, in order to identify about 0.45 km² of unstable areas. Finally, thanks to the availability of ascending and descending results, the time series of deformations of horizontal and vertical components have been obtained using the selected points reported in Figure 6. Within all the UAs identified, six have been selected, based on the highest values of displacement rates, which can be divided into three categories (reported, unknown tunnels and landslides) as shown in Table 2. In order to reconstruct the kinematic of the deformation, vertical and horizontal components have been calculated (

Table 3. Line of sights (LoS), vertical and horizontal displacement components for ascending and descending orbits. (Incident angle: 34° ascending; 43° descending).

	DLoS Ascending (cm)	DLoS Descending (cm)	Dvertical (cm)	Dhorizontal (cm)
UA1	5.0	3.0	4.0	2.0
UA2	2.5	2.0	2.0	0.2
UA3	3.0	2.0	2.5	0.5
UA4	3.5	6.0	6.0	3.0
UA5	2.0	2.5	2.5	0.2
UA6	1.5	2.5	2.5	1.5

) by means of the composition of ascending and descending data for the time span when both acquisitions were present (from January 2016 to September 2018) (Figure 8).

Specifically, the cumulated LoS displacements for the A1/D1, A2/D2 and A3/D3 are about 5.0/−3.0, 2.5/2.0 and 2.5/2.0 cm. The composition of the ascending and descending data made it possible to evaluate the horizontal and vertical components, highlighting in all three cases how the vertical component is always the predominant one, varying between 2.0 and 4.0 cm. This result, combined with the morphology of the sites (flat areas), allows us to state that a subsidence phenomenon is taking place in the area. Moreover, as far as UA2 concerns where the “Inmaculada Fe y Alegria” school was located, it is possible to highlight that the subsidence was already active before the event occurred in 2017 and nowadays it seems to move southward [48]. In the second categories, UA4 and UA5, cumulated LoS displacements are about 3.5 and 2.0 cm in ascending and 6.0 and 2.5 cm for the descending orbit, respectively. Therefore, as previously done, the horizontal and vertical components have been calculated (Table 3). In this case, the vertical components recorded higher values than the horizontal ones, varying between 2.5 and 6.0 cm, but unlike the first three UAs, it can be assumed that rotational landslides are taking place, due to the morphology of the sites (steep slopes). These results are also confirmed by ARCOM [49], which reported an increase in damage to structures in the area after the 2017 earthquake.

Finally, in the last category (UA6), displacements along the LoS of about 1.5 cm in ascending and 2.5 cm in descending, respectively, were recorded. In this case, the vertical component was higher than the horizontal one (Table 3), 2.5 with respect to 1.5 cm, confirming the subsidence phenomenon of the area, even in the absence of reported tunnels.

6. Conclusions

The Zaruma town has been known for ages for gold mineralization and for its mining activity which were considered the main sources of the country economy from the beginning of the 20th century. The activity was divided between small mining concessions released to private local miners, as well as to a few foreign companies. However, simultaneously with the licensed mining activity, a widespread illegal gold exploitation also developed in areas outside the allowed concessions. The forbidden mining area created in 1992, covering 0.7 Km² of the Zaruma urban area, recently has been expanded to 1.73 Km² due to the development of illegal mining activity.

Nowadays, numerous subsidence events occurred in the city, sometimes caused by the mining activities, as the sinkhole generated near the “La Immaculada de fe y alegría” school. Ecuador’s Agency for Mining Control and Regulation (ARCOM) detected in the country about 65 km of tunnel mining, where an unquantified number of illegal miners work, and open, new mines beneath from about 10 to 20 m [50] have been reported. Considering the extension of the underground tunnels, it is difficult to carry out internal controls. According to the Undersecretary of Zaruma Risk Management, the ground surface vulnerability is probably due either to the illegal underground mining or to incorrect urban planning [51].

Thus, it is essential not to underestimate the potential effects of ground surface displacements and it is useful to monitor these to minimize negative environmental and social impacts. Several studies affirmed that in order to decrease the damage due to anthropogenic subsidence phenomena, it is important to identify and to map them [52,53,54,55]. Recently, the management of terrain deformation has increased in importance and several monitoring techniques have been implemented, among which is the differential interferometric synthetic aperture radar (DInSAR) technique.

The DInSAR technique has proved effective in identifying areas subject to subsidence in the Zaruma town, making it possible to recognize critical zones both in areas where the presence of tunnels is known, but especially in those areas where tunnels have never been identified. The latter, which represent a serious problem for the city of Zaruma, as highlighted by several newspaper articles and as the intervention of the Ecuadorian authorities testified, was the main goal and novelty of this work. This result will certainly be very useful to forecast the occurrence of events such as the one recorded in 2017 in the “La Immaculada de fe y alegría” school and in order to help the authorities that have to manage the territory through the identification of the most appropriate risk mitigation actions.

Finally, this result has a double value: on one hand, it will be possible to keep the ground deformation rates under control, in correspondence to the already known tunnels, by carrying out continuous monitoring and integrating it with in situ measurements, thanks to the availability of images of the SENTINEL-1 constellation which have a revisiting time of 12 days; on the other hand, it will be possible to carry out a validation through field detection by focusing attention directly on the areas where anomalous displacement rates have been identified, but where, up to now, there is no certainty of the presence of reported tunnels. Updating the tunnel map, by means of the integration of remote sensing and fields surveys, will therefore be fundamental to subsequently carry out a sinkhole hazard study [56,57,58,59,60], which will allow the authorities in charge of land management to implement risk mitigation actions.