Next Article in Journal
Comparison of Object Detection and Patch-Based Classification Deep Learning Models on Mid- to Late-Season Weed Detection in UAV Imagery
Next Article in Special Issue
Earth Environmental Monitoring Using Multi-Temporal Synthetic Aperture Radar: A Critical Review of Selected Applications
Previous Article in Journal
Attenuation Correction of X-Band Radar Reflectivity Using Adjacent Multiple Microwave Links
Previous Article in Special Issue
A Novel Active Contours Model for Environmental Change Detection from Multitemporal Synthetic Aperture Radar Images

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

Department of Earth Sciences, Environment and Resources, University of Naples, Federico II, 80126 Naples, Italy
Department of Earth Sciences, Natural History Museum, London SW7 5BD, UK
Instituto Ecuatoriano de Regimen Seccional (IERSE), University of Azuay, EC-01.01.981, Cuenca 010107, Ecuador
Author to whom correspondence should be addressed.
Remote Sens. 2020, 12(13), 2135;
Received: 25 May 2020 / Revised: 30 June 2020 / Accepted: 1 July 2020 / Published: 3 July 2020
(This article belongs to the Special Issue Multi-temporal Synthetic Aperture Radar)


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.
Keywords: mining subsidence; DinSAR; risk management mining subsidence; DinSAR; risk management

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), 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:
KDE = 1 ( B a n d w i d t h ) 2 i = 1 n [ 3 π ( w i ) ( 1 ( L i B a n d w i d t h ) 2 ) 2 ]
f o r   L i < B a n d w i d t h
  • Bandwidth is the radius;
  • wi is the weight;
  • Li 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:
B a n d w i d t h = 0.9 × min ( SD , 1 ln ( 2 ) × D m ) × n 0.2
* SD = i = 1 n w i ( x i X w ) 2 i = 1 n w i + i = 1 n w i ( y i Y w ) 2 i = 1 n w i
  • Min is the minimum value between the two in the interlude;
  • SD is the standard distance;
  • Dm is the (weighted) median distance from (weighted) the mean center;
  • n is the number of points;
  • xi, yi are the coordinates;
  • wi 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 (wi) 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.
{ D a = D x s x a s c × D y s y a s c × D z s z a s c D d = D x s x d e s c × D y s y d e s c   × D z s z d e s c
where Da and Dd are, respectively, the displacement values in ascending and descending geometry, Dh and Dv are the displacement vector components along the horizontal (E–W) and vertical directions and Sxasc, Syasc, Szasc, Sxdesc, Sydesc, and Szdesc 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: (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 km2 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) 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 Km2 of the Zaruma urban area, recently has been expanded to 1.73 Km2 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.

Author Contributions

Conceptualization, D.D.M., N.M.; methodology, L.A.; software, L.A. and D.D.M.; validation, L.A., D.D.M.; resources, R.A.R., C.S.; writing—original draft preparation, L.A., R.A.R. and C.S.; writing—review and editing, N.M.; supervision, D.D.M. All authors have read and agreed to the published version of the manuscript.


This research received no external funding.


Many thanks are due to DARES technology for the support to DInSAR elaborations.C.U.G.RI (Consortium between the Federico II University of Naples and the University of Salerno for the Prediction and Prevention of Major Hazards) is gratefully acknowledged for supporting this research with hardware/software facilities. The authors are also grateful to the precious support given by the anonymous reviewers for their valuable observations and suggestions to improve the quality of the paper.

Conflicts of Interest

The authors declare no conflict of interest.


  1. Artisanal and Small-Scale Gold Mining. Available online: (accessed on 11 May 2020).
  2. Artisanal Mining: The Dangers that Come with the Job. Available online: (accessed on 11 May 2020).
  3. Bell, F.; Stacey, T.; Genske, D. Mining subsidence and its effect on the environment: Some differing examples. Environ. Geol. 2000, 40, 135–152. [Google Scholar] [CrossRef]
  4. Gee, D.; Bateson, L.; Sowter, A.; Grebby, S.; Novellino, A.; Cigna, F.; Marsh, S.; Banton, C.; Wyatt, L. Ground motion in areas of abandoned mining: Application of the intermittent SBAS (ISBAS) to the Northumberland and Durham Coalfield, UK. Geosciences 2017, 7, 85. [Google Scholar] [CrossRef]
  5. Gutiérrez, F.; Parise, M.; De Waele, J.; Jourde, H. A review on natural and human-induced geohazards and impacts in karst. Earth-Sci. Rev. 2014, 138, 61–88. [Google Scholar] [CrossRef]
  6. Parise, M. A procedure for evaluating the susceptibility to natural and anthropogenic sinkholes. Georisk Assess. Manag. Risk Eng. Syst. Geohazards 2015, 9, 272–285. [Google Scholar] [CrossRef]
  7. Fazio, N.L.; Perrotti, M.; Lollino, P.; Parise, M.; Vattano, M.; Madonia, G.; Di Maggio, C. A three-dimensional back-analysis of the collapse of an underground cavity in soft rocks. Eng. Geol. 2017, 228, 301–311. [Google Scholar] [CrossRef]
  8. Fiaschi, S.; Tessitore, S.; Bonì, R.; Di Martire, D.; Achilli, V.; Borgstrom, S.; Ibrahim, A.; Floris, M.; Meisina, C.; Ramondini, M.; et al. From ERS-1/2 to Sentinel-1: Two decades of subsidence monitored through A-DInSAR techniques in the Ravenna area (Italy). GISci. Remote Sens. 2017, 54, 305–328. [Google Scholar] [CrossRef]
  9. Fan, H.; Gao, X.; Yang, J.; Deng, K.; Yu, Y. Monitoring mining subsidence using a combination of phase-stacking and offset-tracking methods. Remote Sens. 2015, 7, 9166–9183. [Google Scholar] [CrossRef]
  10. Przyłucka, M.; Herrera, G.; Graniczny, M.; Colombo, D.; Béjar-Pizarro, M. Combination of conventional and advanced DInSAR to monitor very fast mining subsidence with TerraSAR-X Data: Bytom City (Poland). Remote Sens. 2015, 7, 5300–5328. [Google Scholar] [CrossRef]
  11. Du, Z.; Ge, L.; Li, X.; Ng, A.H.M. Subsidence monitoring over the Southern Coalfield, Australia using both L-Band and C-Band SAR time series analysis. Remote Sens. 2016, 8, 543. [Google Scholar] [CrossRef]
  12. Pawluszek-Filipiak, K.; Borkowski, A. Monitoring mining-induced subsidence by integrating differential radar interferometry and persistent scatterer techniques. Eur. J. Remote Sens. 2020, 53, 1–3. [Google Scholar] [CrossRef]
  13. Pawluszek-Filipiak, K.; Borkowski, A. Integration of DInSAR and SBAS Techniques to Determine Mining-Related Deformations Using Sentinel-1 Data: The Case Study of Rydułtowy Mine in Poland. Remote Sens. 2020, 12, 242. [Google Scholar] [CrossRef]
  14. Billingsley, P. Geology of the Zaruma gold district of Ecuador. Am. Inst. Min. Metall. Eng. 1926, 74, 255–275. [Google Scholar]
  15. Cando Jácome, M.; Martinez-Graña, A.M.; Valdés, V. Detection of Terrain Deformations Using InSAR Techniques in Relation to Results on Terrain Subsidence (Ciudad de Zaruma, Ecuador). Remote Sens. 2020, 12, 1598. [Google Scholar] [CrossRef]
  16. Pilatasig, L.; Gordón, D.; Palacios, O.; Sánchez, J. Dominios Litotectónicos Del Sur Del Ecuador y Norte Del Perú; Technical Report; Proyecto Multinacional Andino Geociencias para las Comunidades Andinas-Ecuador-Perú-Canadá, 2015.
  17. Litherland, M. The Metamorphic Belts of Ecuador; British Geological Survey: Nottingham, UK, 1994. [Google Scholar]
  18. Pratt, W.; Figueroa, J.; Flores, B. Geology of the Cordillera Occidental of Ecuador between 3°00′ and 4°00′ S. World Bank Min. Dev. Environ. Control Proj. CODIGEM Rep. 1997, 1, 58. [Google Scholar]
  19. Vikentyev, I.; Banda, R.; Tsepin, A.; Prokofiev, V.; Vikentyeva, O. Mineralogy and formation conditions of Portovelo-Zaruma gold-sulphide vein deposit, Ecuador. Geochem. Mineral. Petrol. 2005, 43, 148–154. [Google Scholar]
  20. Baldock, J. Geology of Ecuador: Explanatory Bulletin of the National Geological Map of the Republic of Ecuador 1:100,000 Scale 1982; Technical Report; Ministerio de Recursos Naturales y Energéticos: San Francisco, Costa Rica, 1982.
  21. Van Thournout, F.; Salemink, J.; Valenzuela, G.; Merlyn, M.; Boven, A.; Muchez, P. Portovelo: A volcanic-hosted epithermal vein-system in Ecuador, South America. Miner. Dep. 1996, 31, 269–276. [Google Scholar] [CrossRef]
  22. Chiaradia, M.; Fontboté, L.; Beate, B. Cenozoic continental arc magmatism and associated mineralization in Ecuador. Miner. Dep. 2004, 39, 204–222. [Google Scholar] [CrossRef]
  23. Hedenquist, J.W.; Arribas, A.R.; Gonzalez-Urien, E. Exploration for epithermal gold deposits. Rev. Econ. Geol. 2000, 13, 45–77. [Google Scholar]
  24. González, A.O.O.; Pozo, A.F.R.; Amaya, R.J.G. Inestabilidad del terreno en zonas de actividad minera: Caso ciudad de Zaruma Ecuador. Redes De Ing. 2017, 8, 69–81. [Google Scholar] [CrossRef]
  25. Murillo, R. ZARUMA, Historia Minera. Identidad en Portovelo; ABYA-YALA: Quito, Ecuador, 2000; Volume 1. [Google Scholar]
  26. GAD ZARUMA on Twitter: “Colapsó Aula afecta por Hundimiento en la Escuela La Inmaculada. Available online: (accessed on 13 May 2020).
  27. URGENTE: Nuevo hundimiento en Zaruma activa al COE cantonal—MachalaMovil. Available online: (accessed on 13 May 2020).
  28. Zaruma se hunde por la minería clandestina. El Comercio. Available online: (accessed on 13 May 2020).
  29. Gobierno Nacional Presentará el Plan Zaruma—Portovelo 2022. Available online: (accessed on 11 May 2020).
  30. Ludizaca Piedra, E.F.; Ochoa Briones, K.H. Propuesta Para Implementación de un Sistema Integrado de Gestión de Riesgos Por Hundimiento de Suelo en la Zona Urbana del Cantón Zaruma; Trabajo de titulación; Environmental engineering, Universidad Politécnica Salesiana: Cuenca, Ecuador, 2018. [Google Scholar]
  31. Franceschetti, G.; Migliaccio, M.; Riccio, D.; Schirinzi, G. SARAS: A synthetic aperture radar (SAR) raw signal simulator. IEEE Trans. Geosci. Remote Sens. 1992, 30, 110–123. [Google Scholar] [CrossRef]
  32. Cigna, F.; Confuorto, P.; Novellino, A.; Tapete, D.; Di Martire, D.; Ramondini, M.; Calcaterra, D.; Plank, S.; Ietto, F.; Brigante, A.; et al. 25 years of satellite InSAR monitoring of ground instability and coastal geohazards in the archaeological site of Capo Colonna, Italy. In Proceedings of the SPIE 10003, SAR Image Analysis, Modeling, and Techniques XVI, Edinburgh, UK, 18 October 2016. [Google Scholar] [CrossRef]
  33. Pappalardo, G.; Mineo, S.; Angrisani, A.C.; Di Martire, D.; Calcaterra, D. Combining field data with infrared thermography and DInSAR surveys to evaluate the activity of landslides: The case study of Randazzo Landslide (NE Sicily). Landslides 2018, 15, 1–21. [Google Scholar] [CrossRef]
  34. Pastor, J.L.; Tomás, R.; Lettieri, L.; Riquelme, A.; Cano, M.; Infante, D.; Ramondini, M.; Di Martire, D. Multi-Source Data Integration to Investigate a Deep-Seated Landslide Affecting a Bridge. Remote Sens. 2019, 11, 1878. [Google Scholar] [CrossRef]
  35. Pepe, S.; De Siena, L.; Barone, A.; Castaldo, R.; D’Auria, L.; Manzo, M.; Casu, M.; Fedi, M.; Lanari, R.; Bianco, F.; et al. Volcanic structures investigation through SAR and seismic interferometric methods: The 2011–2013 Campi Flegrei unrest episode. Remote Sens. Environ. 2019, 234, 111440. [Google Scholar] [CrossRef]
  36. Ferretti, A.; Prati, C.; Rocca, F. Permanent scatterers in SAR interferometry. IEEE Trans. Geosci. Remote Sens. 2001, 39, 8–20. [Google Scholar] [CrossRef]
  37. European Space Agency—ESA, Sentinels Scientific Data Hub. Available online: (accessed on 13 May 2020).
  38. Mora, O.; Mallorqui, J.J.; Broquetas, A. Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images. Geosci. Remote Sens. 2003, 41, 2243–2253. [Google Scholar] [CrossRef]
  39. Iglesias, R.; Mallorqui, J.J.; Monells, D.; López-Martínez, C.; Fabregas, X.; Aguasca, A.; Gili, J.A.; Corominas, J. PSI deformation map retrieval by means of temporal sublook coherence on reduced sets of SAR images. Remote Sens. 2015, 7, 530–563. [Google Scholar] [CrossRef]
  40. Silverman, B.W. Density Estimation for Statistics and Data Analysis; Chapman and Hall: London, UK, 1986. [Google Scholar]
  41. Lu, P.; Casagli, N.; Catani, F.; Tofani, V. Persistent Scatterers Interferometry Hotspot and Cluster Analysis (PSI-HCA) for detection of extremely slow-moving landslides. Int. J. Remote Sens. 2012, 33, 466–489. [Google Scholar] [CrossRef]
  42. Bonì, R.; Pilla, G.; Meisina, C. Methodology for detection and interpretation of ground motion areas with the A-DInSAR time series analysis. Remote Sens. 2016, 8, 686. [Google Scholar] [CrossRef]
  43. Di Martire, D.; Tessitore, S.; Brancato, D.; Ciminelli, M.G.; Costabile, S.; Costantini, M.; Graziano, G.V.; Minati, F.; Ramondini, M.; Calcaterra, D. Landslide detection integrated system (LaDIS) based on in-situ and satellite SAR interferometry measurements. Catena 2016, 137, 406–421. [Google Scholar] [CrossRef]
  44. Solari, L.; Barra, A.; Herrera, G.; Bianchini, S.; Monserrat, O.; Béjar-Pizarro, M.; Crosetto, M.; Sarro, R.; Moretti, S. Fast detection of ground motions on vulnerable elements using Sentinel-1 InSAR data. Geomatics 2018, 9, 152–174. [Google Scholar] [CrossRef]
  45. Solari, L.; Del Soldato, M.; Montalti, R.; Bianchini, S.; Raspini, F.; Thuegaz, P.; Bertolo, D.; Tofani, V.; Casagli, N. A Sentinel-1 based hot-spot analysis: Landslide mapping in north-western Italy. Int. J. Remote Sens. 2019, 40, 7898–7921. [Google Scholar] [CrossRef]
  46. Colesanti, C.; Wasowski, J. Investigating landslides with space-borne Synthetic Aperture Radar (SAR) interferometry. Eng. Geol. 2006, 88, 173–199. [Google Scholar] [CrossRef]
  47. Cascini, L.; Fornaro, G.; Peduto, D. Advanced low-and full-resolution DInSAR map generation for slow-moving landslide analysis at different scales. Eng. Geol. 2010, 112, 29–42. [Google Scholar] [CrossRef]
  48. Morejón, M.; Andrés, E. Aplicación de modelos teóricos para estimar la subsidencia minera que afecta el área de la Unidad Educativa “La Inmaculada Fe y Alegría” cantón Zaruma; Trabajo de titulación; Carrera de Ingeniería en Minas: Quito, Ecuador, 2019. [Google Scholar]
  49. En Zaruma Piden la Reapertura de su Hospital, Cerrado Hace Cinco Meses. Available online: (accessed on 13 May 2020).
  50. ARCOM Realizó en 5 Meses 150 Colapsos Mineros en la Zona de Exclusión de Zaruma. Available online: (accessed on 11 May 2020).
  51. Explotación Minera y Mala Planificación Urbana son los dos Riesgos de Zaruma. Available online: (accessed on 13 May 2020).
  52. Sunwoo, C.; Song, W.K.; Ryu, D.W. A case study of subsidence over an abandoned underground limestone mine. Geosyst. Eng. 2010, 13, 147–152. [Google Scholar] [CrossRef]
  53. Gutiérrez, F. Sinkhole hazards. In Oxford Research Encyclopedia of Natural Hazard Science; 2016; pp. 1–92. Available online: (accessed on 13 May 2020).
  54. Rosdi, M.A.H.M.; Othman, A.N.; Abdul, M.A.M.Z.Z.; Yusoff, Z.M. Sinkhole susceptibility hazard zones using GIS and analytical hierarchical process (AHP): A case study of Kuala Lumpur and Ampang Jaya. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2017, XLII-4/W5, 145–151. [Google Scholar] [CrossRef]
  55. Taheri, K.; Shahabi, H.; Chapi, K.; Shirzadi, A.; Gutiérrez, F.; Khosravi, K. Sinkhole susceptibility mapping: A comparison between Bayes-based machine learning algorithms. Land Degrad. Dev. 2019, 30, 730–745. [Google Scholar] [CrossRef]
  56. Ardau, F.; Balia, R.; Bianco, M.; De Waele, J. Assessment of cover-collapse sinkholes in SW Sardinia (Italy). Geol. Soc. 2007, 279, 47–57. [Google Scholar] [CrossRef]
  57. Intrieri, E.; Fontanelli, K.; Bardi, F.; Marini, F.; Carlà, T.; Pazzi, V.; Fanti, R. Definition of sinkhole triggers and susceptibility based on hydrogeomorphological analyses. Environ. Earth Sci. 2018, 77, 4. [Google Scholar] [CrossRef]
  58. Todd, A.L.; Ivey-Burden, L. A method of mapping sinkhole susceptibility using a geographic information system: A case study for interstates in the karst counties of Virginia. Available online: (accessed on 13 May 2020).
  59. Pellicani, R.; Spilotro, G.; Gutiérrez, F. Susceptibility mapping of instability related to shallow mining cavities in a built-up environment. Eng. Geol. 2017, 217, 81–88. [Google Scholar] [CrossRef]
  60. Parise, M.; Vennari, C. A chronological catalogue of sinkholes in Italy: The first step toward a real evaluation of the sinkhole hazard. In Sinkholes and the Engineering and Environmental Impacts of Karst. In Proceedings of the 13th Multidisciplinary Conference on Sinkholes and the Engineering and Environmental Impacts of Karst, Carlsbad, NM, USA, 6–10 May 2013; Land, L., Doctor, D.H., Stephenson, J.B., Eds.; National Cave and Karst Research Institute, 2013; pp. 383–392. [Google Scholar]
Figure 1. (a) Geological sketch map wih Cross section (A–A’); (b) the study area with the location of pictures (a, b, c) in Figure 2; (c) Simplified cross-section of Portovelo-Zaruma field (modified from [19]). Reference system: WGS84-UTM-17S.
Figure 1. (a) Geological sketch map wih Cross section (A–A’); (b) the study area with the location of pictures (a, b, c) in Figure 2; (c) Simplified cross-section of Portovelo-Zaruma field (modified from [19]). Reference system: WGS84-UTM-17S.
Remotesensing 12 02135 g001
Figure 2. (a) Sinkhole near the school “Inmaculada Fe y Alegria” [26]; (b) subsidence in Zaruma [27]; (c) tunnel under city patrimonial zone [28].
Figure 2. (a) Sinkhole near the school “Inmaculada Fe y Alegria” [26]; (b) subsidence in Zaruma [27]; (c) tunnel under city patrimonial zone [28].
Remotesensing 12 02135 g002
Figure 3. Mining tunnels in the Zaruma urban area (modified from [30]).
Figure 3. Mining tunnels in the Zaruma urban area (modified from [30]).
Remotesensing 12 02135 g003
Figure 4. Interferograms distribution in (a) ascending; (b) descending.
Figure 4. Interferograms distribution in (a) ascending; (b) descending.
Remotesensing 12 02135 g004
Figure 5. Temporal sublook coherence maps: ascending (a), and (c). Stable coherence scatterers (SCS) maps: (b) ascending, (d) descending.
Figure 5. Temporal sublook coherence maps: ascending (a), and (c). Stable coherence scatterers (SCS) maps: (b) ascending, (d) descending.
Remotesensing 12 02135 g005
Figure 6. Mean displacement rate map: (a) ascending orbit; (b) descending orbit. A1–6 and D1–6, ascending and descending, respectively, are points for which vertical and horizontal components have been calculated (see Discussion section).
Figure 6. Mean displacement rate map: (a) ascending orbit; (b) descending orbit. A1–6 and D1–6, ascending and descending, respectively, are points for which vertical and horizontal components have been calculated (see Discussion section).
Remotesensing 12 02135 g006
Figure 7. Kernel density map: (a) ascending orbits; (b) descending orbits.
Figure 7. Kernel density map: (a) ascending orbits; (b) descending orbits.
Remotesensing 12 02135 g007
Figure 8. The time series of cumulated vertical and horizontal components. Displacement standard deviation, for each differential interferometric synthetic aperture radar (DInSAR) measurement, are also reported (black).
Figure 8. The time series of cumulated vertical and horizontal components. Displacement standard deviation, for each differential interferometric synthetic aperture radar (DInSAR) measurement, are also reported (black).
Remotesensing 12 02135 g008
Table 1. Synthetic aperture radar (SAR) data stacks analyzed in this study.
Table 1. Synthetic aperture radar (SAR) data stacks analyzed in this study.
SatelliteOrbitPeriodNr Scenes
Sentinel1 A/BAscending17/01/2016–12/06201989
Sentinel1 A/BDescending05/062015–23/09201851
Table 2. List of unstable areas (UAs) and probable cause of the observed SAR displacement rates.
Table 2. List of unstable areas (UAs) and probable cause of the observed SAR displacement rates.
Unstable AreaProbable Cause
UA1Reported tunnels
UA2Reported tunnels
UA3Reported tunnels
UA4Unknown tunnels/Landslide
UA5Unknown tunnels/Landslide
UA6Unknown tunnels
UA7Reported and unknown tunnels
UA8Unknown tunnels
UA9Unknown tunnels
UA10Unknown tunnels
Table 3. Line of sights (LoS), vertical and horizontal displacement components for ascending and descending orbits. (Incident angle: 34° ascending; 43° descending).
Table 3. Line of sights (LoS), vertical and horizontal displacement components for ascending and descending orbits. (Incident angle: 34° ascending; 43° descending).
DLoS Ascending
DLoS Descending
Back to TopTop