The Detection of Active Sinkholes by Airborne Differential LiDAR DEMs and InSAR Cloud Computing Tools

InSAR (Interferometric Synthetic Aperture Radar) cloud computing and the subtraction of LiDAR (Light Detection and Ranging) DEMs (Digital Elevation Models) are innovative approaches to detect subsidence in karst areas. InSAR cloud computing allows for analyzing C-band Envisat and Sentinel S1 SAR images through web platforms to produce displacement maps of the Earth’s surface in an easy manner. The subtraction of serial LiDAR DEMs results in the same product but with a different level of accuracy and precision than InSAR maps. Here, we analyze the capability of these products to detect active sinkholes in the mantled evaporite karst of the Ebro Valley (NE Spain). We found that the capability of the displacement maps produced with open access, high-resolution airborne LiDAR DEMs was up to four times higher than InSAR displacement maps generated by the Geohazard Exploitation Platform (GEP). Differential LiDAR maps provide accurate information about the location, active sectors, maximum subsidence rate and growing trend of the most rapid and damaging sinkholes. Unfortunately, artifacts and the subsidence detection limit established at −4 cm/yr entailed important limitations in the precise mapping of the sinkhole edges and the detection of slow-moving sinkholes and small collapses. Although InSAR maps provided by GEP show a worse performance when identifying active sinkholes, in some cases they can serve as a complementary technique to overcome LiDAR limitations in urban areas.


Introduction
The monitoring of ground deformation processes using remote sensing and image processing techniques has experienced an exponential increase in the last decade. In karst terrains, stereoscopic aerial photograph interpretation and the use of contour lines have traditionally been the main approaches used to map sinkholes [1][2][3][4]. Nowadays, the mapping process as well as the accuracy and completeness of the sinkhole inventories can be substantially improved using airborne LiDAR (Light Detection and Ranging) data [5]. An important advantage of this remote-sensing technique includes the possibility of filtering the vegetation to produce accurate and high-resolution bare-ground digital elevation models (DEMs) [6,7]. These models can be used to: (a) generate 3D representations of the terrain (e.g., hillshades, red relief image maps) that facilitate the delineation of the sinkholes [8], (b) extract automatically morphometric parameters [9], and (c) apply automated routines for the extraction of closed depressions [10]. Recent processing algorithms and filters favor the automated detection and morphometric characterization of sinkholes, significantly improving the speed and the efficiency of the sinkhole mapping process using LiDAR datasets [6,[11][12][13][14][15]. Nevertheless, the automatic mapping of sinkholes suffers from From the geological perspective, the studied zones are situated in the central sector of the Ebro Cenozoic basin, with subhorizontally lying salt-bearing evaporites (gypsum, anhydrite, halite and glauberite) of the late Oligocene-Miocene Zaragoza Formation [49][50][51]. The soluble bedrock is mainly covered by gravelly terrace deposits, locally thickened by synsedimentary dissolution-induced subsidence and reaching more than 60 m thick [4]. Multiple lines of evidence, including borehole data [51], paleosubsidence structures affecting both the alluvial cover and the bedrock [4], hydrochemical data [52,53], geomorphic evidence including kilometer-scale enclosed basins [54], and subsidence rates [22] indicate that subsidence in the Ebro and Gállego valleys is largely related to the interstratal dissolution of halite and glauberite beds interbedded within the evaporitic bedrock [4,55]. In the Zuera sinkholes, a 60 m deep borehole penetrated gravel-filled voids at 17 m depth and halite layers starting from a depth of 42 m [56]. In the surroundings of Zaragoza city, the altitudinal correlation between groups of thickened terraces and the distribution of glauberite and halite units indicate that the subsidence phenomenon has been induced by dissolution of those salt layers [4,55]. The analyzed sinkholes are typically hundreds of meters in diameter and are dominated by sagging subsidence (i.e., cover and bedrock sagging sinkholes or cover sagging sinkholes). The active depressions typically display continuous and relatively rapid subsidence related to the high solubility of the evaporite bedrock [22,23,32]. Eventually, these large sinkholes can be locally affected by the occurrence of nested catastrophic collapses with high damaging potential. Most of the subsidence damage in the area is related to old sinkholes recognizable in old aerial photographs and historical topographic maps that were filled by man-made deposits and occupied by human structures [3]. The subsidence rates at some of the most active sinkholes were gathered by low-precision leveling profiles [1], cumulative displacement measured in deformed man-made structures of known age [2,3], retrodeformation analysis of trench logs together with geochronological data and historical Remote Sens. 2021, 13, 3261 4 of 23 information [57][58][59][60][61][62], high-precision leveling profiles [22,23], and multi-temporal terrestrial laser scanner data [23,31,32]. In addition, the irrigation of crop fields provides a substantial extra water input to the alluvial-evaporite aquifer system enhancing karstification and subsidence processes [52,53]. The result of the high risk of sinkhole hazards in combination with poor preventive planning is that numerous buildings have been demolished and many transportation infrastructures need continuous remediation works.

Materials and Methods
The study employs a multidisciplinary set of methods and data including conventional time-lapse aerial stereoscopic photography, field mapping, high-precision leveling, airborne LiDAR, and differential SAR interferometry to assess sinkhole activity.
The sinkhole inventory was generated using: (1) old topographic maps from 1960 to 1974, (2) aerial stereoscopic photographs and orthoimages taken in 1946,1956,1986,1996,2002,2009,2012 and 2018, some downloaded from the Spanish National Geographic Institute (IGN), and (3) a 1 m spatial resolution digital elevation model (DEM) derived from the 2010 and 2016 LiDAR data. This inventory was initially published in [3] and has been continuously updated incorporating new sinkhole occurrences and additional information obtained by field mapping and applying hydrochemical, trenching, geophysical and levelling techniques [4,32,[56][57][58][59][60][61][62][63][64][65]. The inventory was completed by a sinkhole-specific study in the most sinkhole-prone area of Zaragoza city (labelled as a detailed area in Figure 1). Here, several field campaigns were carried out between 2014 and 2016 in order to identify evidence of sinkhole-related active deformation on the ground and in human structures. This detailed information helped us to identify active sinkholes, assess the accuracy of the LiDAR-derived displacement data for active-sinkhole detection, and compare these results with the SBAS and FASTVEL algorithms. Note that the inventory of active sinkholes is a minimum record since the lack of field evidence of active deformation does not necessarily mean that the sinkhole is inactive.
Based on their geometry and the mechanisms involved in the ground subsidence [66], the sinkholes were classified into three main genetic types: (a) circular-shaped, collapse sinkholes often between 2 and 10 m across and occasionally up to 40 m in diameter with scarped or overhanging edges that usually occur in a catastrophic way, (b) circular-shaped, sagging and collapse sinkholes tens of meters across, and (c) large, irregularly shaped compound depressions up to 350 m in diameter and 94,000 m 2 in area that result from the coalescence of several smaller sinkholes. The bottoms of some of these depressions are used to reach the water table, hosting lacustrine or palustrine areas.
Open-source LiDAR data of the Ebro Valley were acquired on 4 January 2010 and 22 February 2016 (2238 days; 6.13 years) by the IGN using Leica ALS-60 and ALS-80 laser scanners operated on an aircraft at an altitude of 3500 m. The data that were freely downloaded from the IGN server were already classified according to the American Society for Photogrammetry and Remote Sensing (ASPRS) standard LiDAR Point Classes [67] using the Axelsson [68,69] and NVDI algorithms. Additional information on LiDAR processing and classification can be found in Lorite-Martínez et al. [70]. Following the recommendations presented by Montealegre et al. [71] and Zhao et al. [72] in their reviews on filtering algorithms from open-source LiDAR point clouds, we applied the progressive TIN (Triangulated Irregular Network) densification algorithm (LAStools), the weighted linear least-squares interpolation-based method (FUSION), and the multiscale curvature classification (MCC) with unsuccessful results. All the filters had difficulties on steep slopes over 30 • mislabeling points as nonground, while points returned from small buildings, dense vegetation areas, light posts and shrub cover along irrigation ditches and cropland boundaries are commonly mislabeled as ground surface. Consequently, a buffer of 3 m on irrigation ditches and buildings was removed from the differential raster due to wrong point cloud classification in order to significantly reduce artifacts. An important implication of this filtering process is that the human structures affected by subsidence were undetected by the analysis. The density of the ground measuring points was at least 0.5 points/m 2 with a mean square error in the vertical and horizontal axes of 20 and 10 cm, respectively. DEMs of a 1 m horizontal resolution were produced with the 2010 and 2016 bare-ground LiDAR data using the ArcGIS tool "Las Dataset to Raster" and the triangulation method. These DEMs were then subtracted on a pixel-by-pixel basis for elevation change detection and for active sinkhole mapping. DEM subtraction on confidently known stable areas showed that the vertical accuracy of the vertical displacement model reaches around 25 cm. Consequently, only active sinkholes with subsidence rates over 4 cm/yr (25 cm in ca. 6 years) and associated with a favorable land cover could be detected.
The GEP (Geohazard Themathic Exploitation Platform) is a web-based platform that offers users the possibility to easily request and download analyses of satellite data to monitor geological hazards. In this case, we have run the web-based SBAS and FASTVEL processing approaches that are available with restricted access on the platform. These approaches have been developed by IREA-CNR and TRE-Altamira, respectively, for monitoring the temporal evolution of surface deformations and generating interferogram stacks from a set of Envisat and Sentinel-1 SAR images. The processing chain finally ends with the creation of LOS mean displacement velocity maps (Multitemporal Analyisis Mode) of the ground surface in cm/year by applying advanced DInSAR (SBAS approach) and PSInSAR (FASTVEL approach) techniques. Additional information about the platform, its services, workflow and processing chain is available in Casu et al. Four surface velocity maps with a spatial resolution of 90 by 90 m using the SBAS Sentinel approach by processing 74 archived SLC_IW images of the Sentinel S1 satellite acquired on ascending orbit from 5 November 2014 to 26 January 2021 (track 103), 51 archived SLC_IW images of the Sentinel S1 satellite acquired on descending orbit from 9 June 2016 to 26 January 2021 (track 8), 100 archived SLC_IW images of the Sentinel S1 satellite acquired on ascending orbit from 13 June 2019 to 20 February 2021 (track 30), and 100 archived SLC_IW images of the Sentinel S1 satellite acquired on descending orbit from 12 June 2019 to 19 February 2021 (track 8).

3.
Four surface velocity maps to cover the study area with a spatial resolution of 40 by 40 m using the FASTVEL Sentinel algorithm by processing 400 archived SLC_IW images of the Sentinel S1 acquired on descending orbit from 12 July 2017 to 21 March 2021 (track 8).
The unstable points in the surface velocity maps were selected establishing an average line of sight (LOS) displacement rate threshold of ±2 and ±5 mm/yr for Envisat and Sentinel-1 SAR images, respectively [76]. Tables S1-S3 (supplementary Materials) show additional information on SAR datasets involved in every processing analysis and the thresholds used to run data processing. The vertical deformation rate for the Zuera sinkhole site and Ebro Valley was determined following the methodology described by Bejar et al. [78] from the Envisat and Sentinel SAR images on the ascending and descending orbits and assuming that the north component is negligible [79,80]. Note that the movement in sinkholes associated to subsidence is vertical rather than eastward. The displacement rates along the satellite line-of-sight (vLOS) from ascending and descending datasets were interpolated using the inverse distance weighted (IDW) method (with power 2, search radius 200 m, and one point minimum). The interpolated ascending (vLOSa) and descending (vLOSd) velocities were then used to calculate the vertical components of the velocity, using the raster calculator tool of the open-source software Quantum GIS.
where Ha, Hd, Ea, and Ed represent the direction cosines of the ascending (a) and descending (d) LOS vector and are estimated from the incidence angles (αa and αd) and the LOS azimuths (γa and γd) in degrees (Tables S1-S3

Active Sinkhole Detection by LiDAR Displacement Maps
The criteria to consider a sinkhole as an active subsidence depression rather than an artifact in the DoD mainly depend on the genetic type of the sinkhole. Collapse sinkholes exhibit a well-defined circular subsidence area that indicates the minimum extent of the sinkholes (Figure 2A). Small collapses may be erroneously interpreted as false positive readings related to wrong point classification and a wide variety of alterations (e.g., anthropogenic ground disturbances, vehicles). Consequently, any small collapse had to be validated as a scarped-edge depression in the 2018 orthoimage (25 cm spatial resolution) or in the field. Unfortunately, most of these small collapse sinkholes are commonly filled immediately after their formation, and therefore they were mostly undetected. Large sagging and collapse sinkholes were easily detected in the LiDAR deformation maps. They showed a concentric pattern of ground deformation with increasing displacement towards their center ( Figure 2B,C). Large compound depressions displayed complex spatial deformation patterns since their activity can focus at specific sectors rather than affecting the totality of the depression. They often showed: (1) a cluster of small collapses with a concentric arrangement or along arcuate lines potentially related to precursory evidence of sudden sinkhole formation, (2) linear subsidence troughs that may link two or more small sinkholes, and (3) subcircular and concentric deformation zones coinciding with the most active sectors ( Figure 2D,E). Those readings from the subtraction of the LiDAR-derived DEMs that did not match any of these subsidence arrangements were considered artifacts.

Zuera Sinkholes Site
The Zuera sinkholes site, located 29 km northeast of Zaragoza city (Figure 1), shows two large subsidence zones to the east of the Gállego River on T12 and T11 terraces [81] with areas of 286,000 and 143,000 m 2 , respectively ( Figure 3).
The western sinkhole hosts an ephemeral saline wetland with salt crusts and rudimentary evaporation ponds used in the past for salt mining. Two water samples collected on 13 November 2016 yielded high concentrations of sodium chloride and calcium sulfate with total dissolved solids of 92 and 95 g/L (Table S4, Supplementary Materials). The water samples were saturated with respect to calcite and dolomite, close to saturation with respect to gypsum and anhydrite and still far from the equilibrium with respect to halite despite their high salinity, as determined by the PHREEQC model and Pitzer database [82]. The sodium-chloride hydrochemical facies of the water, unsaturated with respect to halite (saturation index of −1.5), and within a fluvial environment, clearly indicates that it is related to groundwater discharge from the evaporitic aquifer and that salt dissolution is occurring beneath the sinkhole. The subsidence depression shows a N-S elongated shape Remote Sens. 2021, 13, 3261 7 of 23 with an irregular northern margin due to sinkhole coalescence. It hosts several nested collapse sinkholes that contain permanent ponds. The eastern sinkhole is traversed by the A23 highway and the N230 road, which require frequent remediation works. The colonization village of Puilatos was built in 1956 around 500 m north of this sinkhole, was abandoned in 1975, and subsequently demolished due to severe subsidence damage [56]. The sinkhole can be recognized in aerial photographs taken in 1957 as a dry NW-SE-oriented enclosed basin used for cultivation. However, nowadays, because of the rapid subsidence rate, the topographic surface has intercepted the water table of the alluvial aquifer and hosts a permanent lake around 250 m in diameter. A water sample collected on 13 November 2016 yielded a salinity of almost 2 g/L, indicating some provenance from the halite bearing evaporitic aquifer affected by active dissolution (Table S4,

Zuera Sinkholes Site
The Zuera sinkholes site, located 29 km northeast of Zaragoza city (Figure 1), shows two large subsidence zones to the east of the Gállego River on T12 and T11 terraces [81] with areas of 286,000 and 143,000 m 2 , respectively ( Figure 3).

Zuera Sinkholes Site
The Zuera sinkholes site, located 29 km northeast of Zaragoza city (Figure 1), shows two large subsidence zones to the east of the Gállego River on T12 and T11 terraces [81] with areas of 286,000 and 143,000 m 2 , respectively ( Figure 3).  The DoD shows a large number of artifacts associated with cropland boundaries and human excavations (Figure 3). Despite this limitation, it captures the areas affected by the dissolution-induced subsidence that are coherent with the surface deformation features recognized in the field and historical aerial photographs. In the western sinkhole, ground subsidence is detected along a NE-SW belt and in an area associated with the SE sector of the depression. Most of the recorded subsidence rates vary from −4 to −12 cm/yr, although local areas around the nested collapse sinkholes reach rates of up to −20 cm/yr. The eastern sinkhole shows increasing inward ground deformation values reaching −25 cm/yr at the most subsiding sector. The deformation extends across the A23 highway and the N230 road that subside at rates of up to −16 cm/yr, according to high-precision leveling data. Figure 3 shows a LiDAR subsidence rate profile drawn along the N230 road coinciding with the high-precision leveling line installed by Desir et al. [22]. Both vertical rate profiles show a similar pattern and comparable values, except for the most subsiding zone with a reverse shape in the two troughs. Considering the highest vertical displacement along the line and a time lapse of 2236 days between both flights, the maximum subsidence rate can be estimated at around −13 cm/yr, which is consistent with the rate of −16 cm/yr obtained by high-precision leveling [22]. This difference can be attributed to factors such as: (1) the much lower vertical accuracy of the LiDAR data, (2) the lower spatial resolution of the leveling data (10 m of spacing), and (3) the different monitoring periods and spans of both techniques (2236 days for Lidar and 566 days for leveling). Unexpectedly, the deformation pattern extends beyond the mapped southeast sinkhole edge, indicating a larger size for the depression. This ground-lowering trend was undetected in aerial photographs and field surveys and shows the potential effectiveness of the LiDAR data to improve the delineation of the edges of some active sinkholes.
The InSAR velocity maps provided different results. The SBAS Envisat map detected subsidence only on the eastern sector of the western sinkhole with negative vertical velocities lower than −1 cm/yr, while the eastern and fast-moving sinkhole was undetected. SBAS Sentinel maps in ascending and descending geometries and FASTVEL Sentinel maps in the descending geometry reveal maximum negative vertical velocities of −1 to −2 cm/yr on the central sector of the western sinkhole and −2 to −4 cm/yr along the A23 highway and N330 road at the eastern sinkhole. All of the InSAR processing approaches show a loss of coherence in agricultural lands and underestimate the rate and spatial distribution of the subsidence phenomenon. A detailed pixel-by-pixel subsidence rate comparison between the SBAS Sentinel vertical deformation map derived from the combination of Remote Sens. 2021, 13, 3261 9 of 23 the ascending and descending orbits (see Section 3) and the DoD map in a surface area of 9000 m 2 along the A23 highway in the western sector of the sinkhole evidences that the average settlement of the highway was around 7 and 3.5 cm/yr·m 2 for the DoD and InSAR, respectively ( Figure 3). This means that the best InSAR cloud computing performance in the ascending and descending geometries underestimated subsidence by half (Table S5 and Figure S1, Supplementary Materials).

Active Sinkholes in the Ebro River Valley
The study area covers a stretch of 20 km of the Ebro River valley from Zaragoza city to the Jalón River ( Figure 1). In this sector, the Ebro River has excavated a broad alluvial valley in the evaporitic Zaragoza Formation with a markedly asymmetric geometry, characterized by a prominent and linear gypsum escarpment on the northeastern margin and a stepped sequence of terraces and pediments on the opposite side [83]. Sinkhole development mainly occurs on the alluvium-covered areas and especially on the lowermost terraces affected by intensive irrigation [52,53]. The Imperial Canal that provides the majority of the diluted water to the alluvial-karst aquifer system, roughly defines the southern limit of the most sinkhole-prone areas ( Figure 1). The sinkhole inventory map includes a total of 349 sinkholes that represent 1.77% (2.41 km 2 ) of study area.
The performance of the LiDAR-derived displacement data and the InSAR velocity maps for the detection of active sinkholes is illustrated in Figure 1 and Table 1. The DoD suggests that 27.5% of the sinkholes (101 out of 349) show deformation rates over the detection limit established at −4 cm/yr. In contrast, the InSAR velocity maps indicate that a limited number of measuring points consistently indicate unstable areas with negative mean vertical displacement rates. SBAS Envisat, SBAS Sentinel and FASTVEL Sentinel maps permit the detection of 2, 6 and 6.5% of the sinkholes, respectively, with an increasing detection rate as both the time lapse between scenes and the pixel size decrease. Consequently, the FASTVEL processing approach has yielded the best performance, although with slight differences with the SBAS Sentinel approach. The effectiveness of open-source LiDAR data for sinkhole detection was up to 4 times higher than the best InSAR approach (Figure 1). The size of the sinkholes is a critical factor for InSAR detection since only those with unstable surface areas over 2300 m 2 have been identified. This minimum detectable sinkhole area excludes a large proportion of the collapse or sagging sinkholes in the study area (171 out of 349; 48%).   The analysis of the detailed area framed in Figure 1 and depicted in Figure 4 illustrates the possibilities and limitations of the available LiDAR and InSAR data for identifying active sinkholes and assessing their kinematics in agricultural settings and infrastructures. This representative zone, with an area of 23.5 km 2 , is characterized by the occurrence of pre-existing sinkholes buried by man-made deposits that are responsible for the most severe subsidence damage on human structures. Based on conventional mapping using stereoscopic aerial images and fieldwork, 38 sinkholes out of 145 show clear signs of activity such as concentric cracking, damage on human structures or small collapses [3,64]. Obviously, the actual number of active sinkholes can be significantly larger since the conditions required for the identification of instability evidence are not always met.
the most severe subsidence damage on human structures. Based on conventional mapping using stereoscopic aerial images and fieldwork, 38 sinkholes out of 145 show clear signs of activity such as concentric cracking, damage on human structures or small collapses [3,64]. Obviously, the actual number of active sinkholes can be significantly larger since the conditions required for the identification of instability evidence are not always met.  The filtered dense point clouds of the LiDAR data unfortunately show a large number of artifacts that may be misinterpreted as small collapses and require a time-spending validation. However, it successfully identifies most of the active, fast-moving sagging sinkholes and compound depressions. Overall, it allows for the identification of 21 out of the 38 previously mapped active sinkholes of any type (sagging, collapse, or a combination of both) and size (small collapses to very large depressions), which represents a success rate of 57% (Figure 4). The undetected ones may correspond to slow-moving sinkholes with subsidence rates under the threshold limit (−4 cm/yr) and/or sinkholes that have experienced artificial changes between the first and second LiDAR flight.
On the contrary, LiDAR data help us to change the category of 14 sinkholes previously classified as potentially inactive and recognize 5 possible new subsidence areas that were not included in the sinkhole inventory. Subsequent field checking revealed that two out of the five potential sinkholes correspond to true dissolution-induced subsidence depressions, one remains uncertain due to lack of accessibility, and the remaining two are plots of land subject to elevation changes related to human alterations. The addition of detected active sinkholes by conventional mapping plus differential LiDAR equals 54 out of 145 sinkholes (37% of the sinkhole inventory). If we consider the total number of active sinkholes (conventional mapping plus differential LiDAR), the success rate of the DoD reaches 68%. These values reveal that: (1) LiDAR allows for the identification of sinkholes with subsidence rates over −4 cm/yr, (b) the combination of conventional mapping and LiDAR data allows one to expand the number of sinkholes classifiable as active, and (c) almost half of the sinkhole population shows signs of ongoing subsidence.
The InSAR maps generated with the SBAS Envisat, SBAS Sentinel and FASTVEL Sentinel algorithms only provide subsidence data for 6, 11 and 17 active sinkholes that represent 4, 7 and 12% of all the sinkhole inventory, respectively. Considering the totality of active sinkholes (conventional mapping plus LiDAR), the higher spatial resolution of FASTVEL maps showed the better performance, with a success rate of around 30%, although it was less efficient than LiDAR and conventional mapping alone. The lower performance can be attributed to the low InSAR data point coverage in cultivated lands, the small size of many sinkholes and the low spatial resolution of the InSAR maps [64]. In fact, all the sinkholes detected by InSAR correspond to very large compound sinkholes located in urbanized areas or along infrastructures rather than in agricultural terrains, due to the signal decorrelation in vegetated and non-developed environments. The Papiro sinkhole ( Figure 5) is firstly recognized in the 1956 aerial image as a large, 130,000 m 2 compound sinkhole comprising a large central depression and two other sinkholes that hosted several lakes. Since then, this compound sinkhole has been partially filled and several buildings have been built on top. The N232 road that runs across its NE sector has been interrupted several times because of the development of sudden collapses [22]. High-precision leveling data acquired along a 100 m long line with a spacing between benchmarks of 3 m (March 2015 to September 2016) yielded a maximum subsidence rate of −1.9 cm/yr [22]. LiDAR DEM subtraction evidences that the southwestern sinkhole and the large central compound depression are still active with subsidence rates of −6 and more than −25 cm/yr, respectively. The western sector of the larger depression is the most affected by subsidence that shows a distinct linear NW-SE deformation band. Although field evidence, precise leveling data and the development of recurrent collapse sinkholes demonstrate that the N232 road is affected by subsidence, the LiDAR displacement map does not capture deformation. A subsidence-rate profile along the leveling line surveyed by Desir et al. [22] shows a serrated pattern, with apparent uplift in the sinkhole margins and subsidence in the sinkhole with similar rates to those measured by leveling. The deviations can be ascribed to the re-asphalting of the road carried out after the occurrence of a collapse in 2013 ( Figure 5). On the other hand, the deformation map also displays many artifacts mainly located in the northeastern sector of the large depression that are related to densely vegetated areas, a parking lot, cropland boundaries and anthropogenic fillings. Among them, only the two spots located at the northern margin of the compound sinkhole and labelled as 1 and 2 in Figure 5 were validated as small collapse sinkholes of 3 and 5 m in diameter in the 2018 orthoimage, revealing the limitations of open-source, low-resolution airborne LiDAR data for the automatic detection of small collapses. Unexpectedly, none of the InSAR velocity maps identify any ground deformation, pointing to the stability of the sinkhole site and the underestimation of the subsidence hazard. a collapse in 2013 ( Figure 5). On the other hand, the deformation map also displays many artifacts mainly located in the northeastern sector of the large depression that are related to densely vegetated areas, a parking lot, cropland boundaries and anthropogenic fillings. Among them, only the two spots located at the northern margin of the compound sinkhole and labelled as 1 and 2 in Figure 5 were validated as small collapse sinkholes of 3 and 5 m in diameter in the 2018 orthoimage, revealing the limitations of open-source, low-resolution airborne LiDAR data for the automatic detection of small collapses. Unexpectedly, none of the InSAR velocity maps identify any ground deformation, pointing to the stability of the sinkhole site and the underestimation of the subsidence hazard. The Pasarela and Bricodepot sinkholes site includes two very active and damaging sinkholes. The Pasarela sinkhole is a 1000 m 2 , 35 m across, circular depression that is causing severe damage to buildings, a parking lot, and an important water-supply pipe (Figure 6). The Bricodepot sinkhole is a 6500 m 2 , approximately 90 m across, sub-circular depression that affects a department store and its parking lot and the N232 road that requires continuous re-asphalting works. This sinkhole is already recognizable in the 1946 aerial images and was filled between 1957 and 1969 to build an industrial estate and the N232 road on top. Further details on the distribution of surface deformation features can be found in Desir et al. [22] and Benito-Calvo et al. [31]. Soriano and Simón [1] and Simón et al. [2] estimated a maximum subsidence rate of −12 cm/yr for the Bricodepot sinkhole using a home-made leveling device with a bubble level. Subsequently, high-precision leveling data acquired between March 2015 and September 2016 along the northern edge of The Pasarela and Bricodepot sinkholes site includes two very active and damaging sinkholes. The Pasarela sinkhole is a 1000 m 2 , 35 m across, circular depression that is causing severe damage to buildings, a parking lot, and an important water-supply pipe ( Figure 6). The Bricodepot sinkhole is a 6500 m 2 , approximately 90 m across, sub-circular depression that affects a department store and its parking lot and the N232 road that requires continuous re-asphalting works. This sinkhole is already recognizable in the 1946 aerial images and was filled between 1957 and 1969 to build an industrial estate and the N232 road on top. Further details on the distribution of surface deformation features can be found in Desir et al. [22] and Benito-Calvo et al. [31]. Soriano and Simón [1] and Simón et al. [2] estimated a maximum subsidence rate of −12 cm/yr for the Bricodepot sinkhole using a home-made leveling device with a bubble level. Subsequently, high-precision leveling data acquired between March 2015 and September 2016 along the northern edge of the N232 road yielded a reliable maximum subsidence rate of −7.5 cm/yr [22]. Finally, the comparison of 3D point clouds obtained by a terrestrial laser scanner (TLS) every 4 months from October 2014 to April 2017 revealed subsidence with a complex spatial-temporal distribution with maximum subsidence rates of −15 and −8 cm/yr for the Pasarela and Bricodepot sinkholes, respectively, in agreement with the leveling data [31]. the N232 road yielded a reliable maximum subsidence rate of −7.5 cm/yr [22]. Finally, the comparison of 3D point clouds obtained by a terrestrial laser scanner (TLS) every 4 months from October 2014 to April 2017 revealed subsidence with a complex spatial-temporal distribution with maximum subsidence rates of −15 and −8 cm/yr for the Pasarela and Bricodepot sinkholes, respectively, in agreement with the leveling data [31]. The DoD captures surface displacement in both sinkholes. In the Pasarela sinkhole, surface settlement seems to be concentrated in the parking lot and shows a concentric displacement pattern increasing towards the center. In contrast, in the Bricopedot sinkhole, the displacement model indicates subsidence mainly along the NE service road, The DoD captures surface displacement in both sinkholes. In the Pasarela sinkhole, surface settlement seems to be concentrated in the parking lot and shows a concentric displacement pattern increasing towards the center. In contrast, in the Bricopedot sinkhole, the displacement model indicates subsidence mainly along the NE service road, where re-asphalting is carried out with a much lower frequency than that along the main carriageways. The subsidence-rate profiles generated from the displacement models provide maximum subsidence rates of −15 and −7.5 cm/yr for the Pasarela and Bricodepot sinkholes, which are coherent with the already published precise leveling and TLS data [22,31]. The subsidence-rate profile of the Bricodepot sinkhole shows higher subsidence rates along most of the profile than those obtained by high-precision leveling, suggesting a probable systematic error in the elevation data of the filtered LiDAR data. Unfortunately, the LiDAR data only detects ground deformation around the central sector of the sinkholes affected by higher subsidence. However, field mapping, leveling data, and TLS deformation maps indicate that subsidence extends towards the sinkhole edges, although at progressively lower rates. This important limitation is related to the detection limit of −4 cm/yr established for the LiDAR data due to vertical precision limitations. Artifacts are associated with vehicles and dense vegetation and appear as small patches distinguishable from the consistent deformation patterns of the sinkholes ( Figure 6). Once again, the SBAS Envisat, and SBAS and FASTVEL Sentinel processing approaches were not able to identify ground deformation in any of the sinkholes. However, the FASTVEL map depicted three deformation points outside the southwestern boundary of the Bricodepot sinkhole that might be related to incipient subsidence and sinkhole growing direction.

Measuring Subsidence Rates in Active Sinkholes
Finally, the Estrellas sinkhole is a 2000 m 2 , 120 m across, sub-circular depression located in a suburb of Zaragoza city (Figure 7). It has a long history (since the 1960s) of anthropogenic filling, building construction, damage, and demolition [3]. At the present time, a multi-story building at the eastern sector is affected by severe damage [65]. The age of the anthropogenic sinkhole fillings exposed in two trenches excavated in the central sector of the sinkhole across fault scarps results in a minimum average displacement rate between −1.2 and −1.8 cm/yr for a poorly constrained time lapse between 1956 and 1997. Moreover, the vertical displacement of human structures allowed Gutiérrez et al. [57] to estimate rough subsidence rates of −3 cm/year. High-precision leveling data acquired between January 2017 and February 2019 across the sinkhole along an east-west direction yielded a maximum subsidence rate of around −2.5 cm/yr [65]. Neither the DoD nor any of the InSAR velocity maps detected deformation, erroneously suggesting that the area is stable. However, the subsidence-rate profile across this sinkhole yields a maximum lowering rate of around −3.5 cm/yr, which compares well with the leveling data and displacement measurements from human structures, although with a contrasted different geometry and an overestimation of the subsidence in its western sector related to the limitations of the LiDAR-derived displacement models. The Estrellas sinkhole is a good example of the drawbacks of the DoD for the detection of slow-moving active sinkholes in urban areas where the filtering of the point clouds in order to reduce artifacts obliterates relevant subsidence information. Remote Sens. 2021, 13, x FOR PEER REVIEW 15 of 22

Advantages and Limitations of Airborne LiDAR and Web-Based InSAR Processing Tools
The results presented in this paper highlight the potential utility of multi-temporal airborne LiDAR data for the detection of active sinkholes and the assessment of subsidence rates. Nevertheless, there were several limitations derived from the acquisition of the LiDAR data. Firstly, the quality of the LiDAR-derived displacement maps depends on the land cover, offering better results in agricultural terrains and along infrastructures, rather than urbanized areas. Secondly, the low precision, with a vertical displacement threshold of ±25 cm only allows for the detection of fast-moving sinkholes with subsidence rates over −4 cm/yr. In the closest area to Zaragoza city, the proportion of active sinkholes captured by differential LiDAR was 21 out of 38 (success rate of 57%), which represents only half of the previously mapped active sinkholes using conventional methods. On the contrary, it permitted the identification of the most damaging sinkholes, the discovery of 2 new non-inventoried active sinkholes and the re-classification of 14 sinkholes as active. In addition, the effectiveness of multi-temporal LiDAR datasets to roughly estimate subsidence rates is demonstrated by the consistency between the subsidence-rate profiles derived from LiDAR and from high-precision leveling data. It is important to stress that the approach has been applied in a very favorable area in which a significant proportion of the sinkholes are affected by rapid and progressive sagging subsidence related to the dissolution of salt-bearing evaporites. The performance of the method would be much worse in karst areas with slow-moving sinkholes and/or characterized by episodic displacement, which is a more common situation (e.g., carbonate karsts).
The type and size of the sinkholes had a significant impact on the results. Active sagging sinkholes and large compound depressions are identified thanks to their concentric

Advantages and Limitations of Airborne LiDAR and Web-Based InSAR Processing Tools
The results presented in this paper highlight the potential utility of multi-temporal airborne LiDAR data for the detection of active sinkholes and the assessment of subsidence rates. Nevertheless, there were several limitations derived from the acquisition of the LiDAR data. Firstly, the quality of the LiDAR-derived displacement maps depends on the land cover, offering better results in agricultural terrains and along infrastructures, rather than urbanized areas. Secondly, the low precision, with a vertical displacement threshold of ±25 cm only allows for the detection of fast-moving sinkholes with subsidence rates over −4 cm/yr. In the closest area to Zaragoza city, the proportion of active sinkholes captured by differential LiDAR was 21 out of 38 (success rate of 57%), which represents only half of the previously mapped active sinkholes using conventional methods. On the contrary, it permitted the identification of the most damaging sinkholes, the discovery of 2 new noninventoried active sinkholes and the re-classification of 14 sinkholes as active. In addition, the effectiveness of multi-temporal LiDAR datasets to roughly estimate subsidence rates is demonstrated by the consistency between the subsidence-rate profiles derived from LiDAR and from high-precision leveling data. It is important to stress that the approach has been applied in a very favorable area in which a significant proportion of the sinkholes are affected by rapid and progressive sagging subsidence related to the dissolution of salt-bearing evaporites. The performance of the method would be much worse in karst areas with slow-moving sinkholes and/or characterized by episodic displacement, which is a more common situation (e.g., carbonate karsts).
The type and size of the sinkholes had a significant impact on the results. Active sagging sinkholes and large compound depressions are identified thanks to their concen-tric deformation pattern with increasing subsidence rates towards their center and linear troughs. However, the technique displays significant limitations in the identification of small collapses. The large number of artifacts associated to vehicles, densely vegetated areas, steep slopes, road embankments, cropland boundaries and anthropogenic fillings led to numerous false positives and an overestimation of the number of potential collapse sinkholes. Consequently, every small site with apparent subsidence had to be validated with the help of high-precision orthoimages or fieldwork. This is very time-consuming work that requires a specialist with knowledge of the subsidence phenomena and, specifically, the existence of a previous cartographic sinkhole inventory.
The results obtained from the InSAR cloud computing approaches only provide ground motion due to karst subsidence in some specific cases, either in ascending or descending orbits. Among the different algorithms, the FASTVEL maps provided the best detection rate around Zaragoza city (30% of all the active sinkholes) and was the only InSAR approach able to capture displacement within some of the most damaging sinkholes of the analyzed area. These results indicate that the FASTVEL approach yielded the best performance among the GEP web-based InSAR tools, but it was still low compared with conventional sinkhole mapping and differential LiDAR. Considering that subsidence rates in many sinkholes do not vary significantly with time due to their sagging behavior [1,2,22,23,31,65], the success of each InSAR web-algorithm in sinkhole detection is apparently dependent on the processing of the SAR images, rather than on the monitoring period. Based on the characteristics described for every processing approach (Tables S1-S3, supplementary Materials), the processing interferometry technique (DInSAR versus PSIn-SAR) is not a critical variable. In contrast, as expected, sinkhole detection slightly improves with the finer spatial resolution of the SAR images and mainly with the frequency of acquisitions (shorter time interval between images) that allows for reducing the decorrelation between images [84]. Based on our results, large-temporal SAR interferograms cannot measure those sinkholes with the highest subsidence rates. Previous DInSAR velocity maps of the study area with a pixel size between 20 to 90 m obtained from 27 ERS, 29 Envisat and 13 ALOS SAR images processed by the SBAS and SPN (Stable point Network) approaches captured between 4 and 25% of the active sinkholes [58]. These values, which are similar to the ones obtained in this study, highlight the difficulties of InSAR techniques in easily identifying active sinkholes in the mantled karst of the Ebro Valley. Several constrains inherent to InSAR data may explain this very low detection proportion in the study area: (a) Decorrelation in agricultural terrains due to vegetation and surface changes. In fact, the majority of the identified sinkholes were located in urbanized areas. This is an important factor in the study area since irrigation lands cover more than 80% of the Ebro River valley. (b) The low spatial resolution of the SAR models that precludes the detection of the small sinkholes. A 20 m spatial resolution SAR image is unlikely to confidently detect sinkholes less than 100 m in diameter [85,86], and recent studies demonstrate that only SAR images with a very high resolution (pixel resolution ranging between 25 cm and 3 m of TerraSAR-X and COSMO-SkyMed satellite images) are able to identify sinkholes several meters across [25,87,88]. Unfortunately, the 40 m pixel resolution velocity maps of the FASTVEL approach that shows the better performance were useful only for the detection of unstable surface areas over 2300 m 2 . This minimum detectable sinkhole area excludes almost half of the small collapse and sagging sinkholes of the study area (171 out of 349; 48%), and so only the largest ones were detected. (c) A significant loss of coherence over time caused by severe ground deformation related to fast-moving sinkholes [89,90]. Thus, although InSAR cloud computing approaches are a great advance in the analysis of Earth surface processes, our results indicate that these tools are not as effective as other methods in the Ebro Valley, where resources such as LiDAR DEMs and historical aerial photographs and topographic maps are available.

Implications of InSAR Cloud Computing and Differential LiDAR in Sinkhole Risk Management
Subsidence damage in the study area mainly occurs within the boundaries of preexisting sinkholes that depict continuous long-sustained subsidence activity [3,58,59,61]. Therefore, the exclusion for the construction of the areas affected by subsidence and their vicinity would be the most effective mitigation measure. Unfortunately, this preventive planning strategy was not implemented in the study area, and a large number of sinkholes have been filled and obliterated by human activities to build infrastructures and buildings. In these cases, our results support the fact that differential LiDAR on regional and local scales is more advantageous than InSAR cloud computing products when detecting dissolutioninduced subsidence mainly due to the higher spatial resolution of the LiDAR displacement maps. Consequently, differential LiDAR could be a useful tool for sinkhole risk assessment and management. Unlike conventional techniques that require the performance of timeconsuming tasks (e.g., mapping, geodetic surveying), multi-temporal LiDAR datasets can remotely provide sinkhole-related displacement data at a relatively low cost over large areas, such as: (1) the location of the most rapid and damaging sinkholes, (2) the highest subsidence rates of particular sinkholes, (3) the spatial variability of subsidence within the large sinkholes and compound depressions differencing the active and inactive sectors, (4) the potential growing trend of sinkholes outside their boundaries that might forecast future damaged areas, and (5) the identification of new and non-inventoried sinkholes.
Unfortunately, the vertical accuracy inherent to the altitude acquisition of the airborne LiDAR data prevents the detection of slow-moving sinkholes and the precise mapping of the sinkhole edges. In most cases, differential LiDAR provides an incomplete picture of the active subsidence areas. Those subsidence areas affected by lower subsidence rates than the LiDAR detection threshold entail a hazard underestimation. The results indicate that LiDAR displacement models overlook around 40% of the inventoried active sinkholes due to their low ground motion rate or the existence of anthropogenic fillings. Consequently, LiDAR should not be considered an alternative to classical mapping methods, but a highly useful complementary approach to improve and upgrade sinkhole inventories and identify the most-damaging subsidence areas. On the other hand, LiDAR deformation maps compared with TLS, leveling or trenching techniques do not define the boundaries of the areas affected by subsidence, but only those sectors with adequate conditions and moving over the vertical detection limit. These maps clearly underestimate the set-back distance for any particular sinkhole that is a critical parameter for land-use planning [59], and they should be complemented by other more precise ground-based techniques. Nevertheless, these differential LiDAR constraints on sinkhole activity assessment could be substantially overcome through: (1) The increase in point densities by reducing the altitude of the flights, at the expense of a higher cost. This would increase the vertical accuracy of the data and would decrease the detection threshold, favoring the identification of slow-moving sinkholes. (2) The acquisition of flights with a shorter re-visiting time in order to recognize ground surface alterations caused by human activities and assess the potential spatiotemporal changes in sinkhole subsidence rates. (3) The improvement of open-source point cloud classification algorithms that might reduce the amount of artifacts to upgrade the automatic identification of small collapses and the detection of sinkholes in urbanized and densely vegetated areas. Unfortunately, although open-source filtering algorithms have a good overall performance with around an 80% success rate [71,72], their application for sinkhole detection still leads to a large number of false positives and negatives.
Given the limitations of differential LiDAR to detect slow-moving sinkholes and measure subsidence in urban areas where InSAR can have a high performance, the application of both techniques can be compatible and InSAR analyses may complete differential LiDAR studies processing high-resolution X-band SAR images. Since the GEP results almost reach the detection capability of previous InSAR analysis using Envisat images at full resolution [58], InSAR cloud computing can be a useful tool to firstly explore karst regions to get a first impression of sinkhole activity where there is not a deep knowledge of karst subsidence. However, based on the InSAR web-tools performance in the Ebro Valley, the detected subsidence in other areas might only represent a small part of the real situation mainly because of the low spatial resolution of the InSAR velocity maps. Thus, InSAR cloud computing results should not be directly used to assess sinkhole hazards on a regional or local scale alone but they can be a good initial reference for further detailed studies using other methodologies.

Conclusions
The subtraction of DEMs derived from multi-temporal, open-access, low-resolution airborne LiDAR data has resulted in a very useful and cost-effective methodology for the detection of active sinkholes in the mantled evaporite karst of the Ebro Valley (NE Spain). Active sagging sinkholes and large compound depressions are easily identified by differential LiDAR thanks to their concentric deformation pattern with increasing subsidence rates towards their center and/or the occurrence of linear troughs. In these types of sinkholes, LiDAR vertical displacement maps provided important information about the location of the most active sectors, their growing trend and their highest subsidence rates based on the consistency between leveling and LiDAR subsidence profiles. Unfortunately, the technique displays significant limitations in the identification of small collapses that are normally overlooked and often interpreted as artifacts, the detection of slow-moving sinkholes and the precise mapping of the sinkhole edges. Despite these limitations, the DoD maps show that around 30% of the sinkhole population in a 135 km 2 sector of the Ebro Valley (101 out of 349) shows signs of ongoing subsidence with subsidence rates over −4 cm/yr. These results support the implementation of free-source airborne LiDAR in sinkhole hazard mitigation programs in evaporite karst terrains where a significant proportion of the sinkholes are affected by rapid and progressive sagging subsidence.
SBAS Envisat, SBAS Sentinel and FASTVEL Sentinel maps permit the detection of 2, 6 and 6.5% of the sinkhole population, respectively, with an increasing detection rate, as both the time lapse between scenes and the pixel size decrease. Consequently, the FASTVEL processing approach yielded the best performance, although with slight differences compared to the SBAS Sentinel approach. These results evidence that the capability of airborne LiDAR data for sinkhole detection in the study area was up to four times higher than the best performance of the FASTVEL and SBAS Sentinel approaches that underestimated the number of active sinkholes. In addition, a pixel-by-pixel comparison between the DoD and the SBAS Sentinel vertical deformation maps in the Zuera Sinkholes site shows that the best GEP InSAR performance underestimates subsidence rate by a factor of two. The low detection proportion of InSAR cloud computing tools is mainly related to the low spatial resolution of the SAR models that precludes the detection of the small sinkholes. Thus, although InSAR cloud computing approaches are a great advance in the analysis of Earth surface processes and can complement LiDAR data to detect slight subsidence areas (mm/yr) in urban areas, our results indicate that, as of today, Web-InSAR tools are not as effective as differential LiDAR or conventional mapping for the detection of active sinkholes.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/rs13163261/s1, Figure S1: Vertical deformation rate of the InSAR data obtained from the 2019-2021 sentinel-1 images and the SBAS Sentinel approach Figure S2: LOS velocity and vertical deformation rates of the detailed area around Zaragoza city from the 2003-2010 Envisat SAR images and the SBAS Envisat approach, Figure S3: LOS velocity and vertical deformation rates of the detailed area around Zaragoza city from the 2019-2021 Sentinel-1 SAR images and the FASTVEL Sentinel approach, Table S1: Main characteristics of the SAR Envisat Dataset processing and surface velocity maps, Table S2: Main characteristics of the SAR SBAS Sentinel Dataset processing and surface velocity maps, Table S3: Main characteristics of the SAR FASTVEL Sentinel Dataset processing and surface velocity maps, Table S4: Chemical composition, total dissolved solids (TDS), alkalinity (Alk) and saturation indexes (si) performed with the PHREEQC software and the Pitzer database for the Zuera