1. Introduction
Groundwater resources are very vulnerable due to natural and anthropogenic impacts, especially in karstic terrains of arid and semi-arid zones, due to increased requirements caused by climate change and population expansion [
1,
2]. Karst landscapes are considered highly fragile settings that are severely vulnerable to sinkhole formation. Sinkhole landforms are depression areas with subcircular shapes on a surface typically common in carbonate and/or evaporate rocks; these terrains are characterized by minerals that might dissolve in water and are controlled by (1) texture, (2) mechanical strength (i.e., rocks have lower mechanical strength and have faster and higher solubility behavior), and (3) solubility mechanisms, etc. [
3]. Typically, water molecules absorb atmospheric carbon dioxide (
) and transform it into carbonic acid (
) in soils, which pushes carbonate minerals into the leaching process. The acidic solution
filters through the weak areas (joints, faults, and fractures) and involves the formation of underground cavities [
4,
5].
Land subsidence issues have become increasingly prominent in recent years in many parts of the world. Land subsidence is the lowering of the ground surface elevation due to underground material movement. It occurs due to a complex process controlled by natural or anthropogenic causes or a combination of both effects. Natural causes could be events such as earthquakes, soil compaction, glacial isostatic adjustment, erosion, and sinkhole formation; anthropogenic effects are activities such as the removal of water, oil, natural gas, or mineral resources from the ground by pumping, fracking, or mining activities. Typically, sinkhole collapse happens when the piezometric level goes down. The damage could vary from one region to another and can be affected indirectly by infrastructure presence and population concentrations [
5]. Because of urban expansion and climate change issues, sinkhole subsidence geohazards have been reported worldwide, including in Mexico [
6], the USA [
7], Iran [
8], the UAE [
9], and Turkey [
10].
Subsidence is a common geological hazard and can cause infrastructure damage, ground surface ruptures, increasing flood risk, and adverse socioeconomic impacts on communities. Detecting the spatial extent and monitoring temporal evolution is crucial to determine the causes of subsidence and prevent activities to mitigate the negative effects. Elevation or elevation-change measurements are fundamental to monitoring land subsidence and have been measured by using traditional methods; Global Positioning System (GPS) stations, leveling networks, and recently, satellite InSAR (interferometric synthetic aperture radar) remote sensing are prominent techniques.
In recent years, InSAR has become a powerful remote sensing technique for monitoring global earth changes. Persistent scatterer interferometry synthetic aperture radar (PS-InSAR), which belongs to differential interferometric synthetic aperture radar (DInSAR) groups, is an innovative technology for surface displacement monitoring. These methods are less affected by time and space decoherence and have low-cost advantages and high precision. Phase interferometry contains artifacts that reduce the quality of the interferograms. A small space and temporal baseline of the C-band images were constructed to avoid and minimize errors conserved during the acquisition time of each image [
11]. However, it is the contribution of atmospheric conditions such as temperature, vapor content, and pressure that are converted into small-magnitude errors in interferograms [
11,
12].
The PS-InSAR approach was first realized by Ferretti [
13]. The algorithms identify high persistent scatterer (PS) pixels less affected by time and space decoherence, as well as relatively stable pixels, by analysis of their amplitude scintillations in interferogram series and time/frequency characteristics of each phase. This algorithm works better in urban areas with a large number of artificial man-made structures that increase the identification of scatterers. Persistent scatterer interferometry (PSI) uses phase analysis for PS pixel identification proposed by Hooper [
14]. This algorithm uses a combination of amplitude dispersion index and phases spatial correlation to identify PS points to produce a time series of deformation using the spatially correlated nature, rather than requiring a known temporal dependence [
14].
The present work aims to study the subsidence in a karst semi-arid area that suffers from severe drought and groundwater exploitation problems. Interferometry techniques with large-scale acquisition and high accuracy provide subsidence rates and help to extract and delineate the ground movement changes registered over the study period. The spatiotemporal land deformation distribution identified by exploring InSAR and GPS highlights the serious land subsidence damages that have occurred in arid and semi-arid regions. The findings will help the local authority to improve a strategy for regional deformation and the implementation of inventory maps to avoid human and economic losses.
2. Study Area
The Cheria basin is situated in the southwestern part of Tebessa within the northeast province of Algeria between the latitudes of 35°18′41,62″ and 35°14′51,30″N and longitudes of 7°42′51,40″ and 7°47′34,62”E (
Figure 1). The basin is surrounded by mountain ranges with an elevation exceeding 1000 m, with some areas at elevations of up to 1200 m, such as Djebal Achour and Djebal Dokkan in the east; Djebal Zora and Djebal Boukammech in the south; Djebal Troubia and Djebal Tazbent in the north; and Djebal Kamallel and Djebel El Abtine in the west.
The local geology consists of lithostratigraphic units that vary from bottom to top by Maastrichtian limestone, Eocene limestone, and Mio-Plio-Quaternary deposits (
Figure 2) [
15]. The sinkhole studied is a cover-collapse type characterized by Eocene limestone layers covered with Quaternary deposits (gravel, sand, silt, and clay) (see
Figure 2). They are mainly synsedimentary subsidence types caused by the dissolution of the karstic bedrock [
16]. The tectonic system dissected by conjugation faults is controlled by the local hydrogeographic networks, the NE–SW and the vertical NW–SE diagonal system, and the N–S (Alpine phase) system [
17]. The region has a semiarid climate, and the local precipitation is approximately 340–400 mm/year.
The dominant type of sinkhole present in the study region is the cover-collapse sinkhole. Cover-collapse sinkholes are one of the most dangerous types that can collapse in minutes or even seconds, causing disastrous damages. Real damage caused by past sinkhole collapse events is displayed in
Figure 3.
3. Data and Methods
In this paper, we focus on investigating the subsidence phenomenon and evolution between specific periods selected based on GPS measurements and PS-InSAR techniques. SAR images were obtained from the Sentinel 1-A and 1-B satellites between 2016 and 2022. The SAR scenes selected for the case study were carefully chosen (e.g., based on good weather, etc.). The study area is mostly bare land with a very low percentage of permanent vegetation, which makes it easier to obtain good SAR scenes. The Sentinel-1 short wavelength is sensitive to vegetation, which leads to serious decorrelated noise due to vegetation changes during the revisiting period of the SAR images. The GPS datasets were measured by a Leica 1200 GNSS system receiver with an ATX1230 antenna and an RTX1250GG receiver.
More information about the satellite image dataset and baseline distribution related to acquisition time is provided in
Table 1.
The phase difference between SAR images consists of measuring the ground difference along the radar line of sight (LOS) direction. The interferometric phase comprises the following contribution in Equation (1):
where
represents the residual orbit error caused by an inaccurate orbit sensor;
is related to the topographic variation enhanced with physical features of the ground surface;
is the unwrapping error and the main algorithm process applied to recovering an unambiguous phase and correcting the residual phases (
,
, and
) for accurate high deformation results;
is the residual thermal noise effect; and
is the atmospheric errors during SAR image acquisition that are influenced by various signals. It is the main artifact source and is generated when the microwave passes through the tropospheric layer (i.e., elevation component due to atmospheric vertical stratification, and turbulent component referring to horizontal changes in water vapor distribution over short time intervals) [
11,
18].
Data preparation consisted of the creation of single-look master interferograms from the N + 1 single-look complex (SLC) SAR images with the most optimal configuration. All the other images were denoted as slaves and co-registered to generate interferograms. Integration of external DEM SRTM 1 sec data was completed as a topographic reference, and these data were also used to eliminate the topographic effects and geocoding. All preprocessing steps were carried out using Gamma software (20220701) and MATLAB (9.13 R2022b). The flowchart of the overall methodology is presented in
Figure 4.
The PSI method processed to find pixels not affected by noise with high phase stability is performed using the following steps:
In preliminary analysis loading, persistent scatterer (PS) candidate points are selected as pixels with a value of the amplitude dispersion index (ADI) that is smaller than a threshold.
Estimate phase noise means the atmospheric phase screen (APS) value is contained on each candidate pixel in the interferogram, defined by the spatially correlated phase and uncorrelated terrain errors. For good results, various spatiotemporal filters are used to correct APS and achieve only the deformation part.
Persistent scatterer points are selected according to the atmospheric phase screen (APS) correction parameter and the percentage of random pixels in a scene per density is estimated by application of a probability statistics method.
The PSs selected in the previous step are weeded, removing those that are deemed too noisy due to signal contributions from neighboring ground resolution elements.
The wrapped phase of the selected pixels is corrected for a spatially uncorrelated look angle DEM error.
Three-dimensional unwrapping of the above-mentioned corrected phase PS result is used; unwrapping errors are more likely to occur in a longer perpendicular baseline interferogram.
A spatially uncorrelated look angle SCLA error was calculated in step iii and removed in step v; in step vii, a spatial look angle error is calculated which is due almost exclusively to a spatially correlated DEM error (this includes an error in the DEM itself and incorrect mapping of the DEM into radar coordinates). The master atmosphere and orbit error phase are estimated simultaneously.
Atmospheric filtering and estimation of other spatial correlation error terms are conducted. The results are a data file containing final PS points with a deformation velocity in the precision of mm/year representing the land deformation model of the area of interest [
14].
To compare these measurements with GPS vertical measurements, we transformed SAR measurements recorded along the radar line of sight (LOS) into 3D displacement. The velocity (
) measured can be constructed from SAR image acquisition geometry using Equation (2):
where θ and α are pixel-based radar incident angles and azimuth angles of the satellite, respectively;
,
, and
represent real vertical north–south direction and west–east direction deformations, respectively. If
and
are the velocities along the radar line of sight for ascending and descending tracks, Equation (2) can be rewritten into Equation (3):
To further simplify Equation (3) with three variables, displacements occurring along the north–south direction cannot be measured accurately (almost parallel to the satellite orbit). Assuming that the north–south component projection of the velocity along the LOS radar is negligible for both ascending and descending orbit tracks, in this case, Equation (3) can be approximated and further simplified:
In situ dataset measurements were performed in this region for six years by manually monitoring millimeter surface displacement in this study area (unfortunately, no permanent GPS monitoring station exists). A total of five ground control points (GCP) were selected and used for validation purposes of persistent scatterer results. The GPS observation measurements were obtained with a session length of 30 min for five selected GCP points and 10 min for each of the remaining sites (Profiles AB-CD). The periodic measurements over an extended period are a great approach to flattening the error curve and eliminating outliers and artifacts in GPS measurements [
19] (i.e., the same concept used by InSAR time-series analysis). The Geodetic Reference System 1980 (GRS80) was used for the ellipsoidal heights. Furthermore, the geoid model of EGG2008 was used to convert the ellipsoidal heights into orthometric heights.
It is important to keep in mind that GPS measurement can obtain the Earth’s surface motion information with high precision and can directly reflect the vertical variation characteristics of the observed object. However, measurements cannot be densely organized on a large scale due to the critical limitations of GPS such as low spatial resolution and high cost and time consumption [
20,
21]. Therefore, the PSI results were used and validated with five monitored GPS ground control points received during the period of analysis.
Finally, to assess the discrepancy found between PS-InSAR and GPS, root mean square error (RMSE) was used. RMSE is common and is considered an excellent general-purpose error metric that can be used to (1) indicate the absolute fit of the model to the data and (2) provide the average model prediction error in units of the variable of interest; RMSE is a negatively oriented score, which means lower values are better [
22,
23,
24]. RMSE can be calculated using Equation (5):
where
is the GPS measurements,
is the PS-InSAR measurements, and
is the number of observations available for analysis.
4. Results
For the deformation monitoring in the study area, we used PS-InSAR in combination with GPS measurements, as described in
Section 3. The mean velocity maps of the final geocoded vertical and horizontal displacement generated from Sentinel-1 data are presented in
Figure 5. Color ramp distribution of the velocity maps gives a visual preliminary overview of the dominant movement (e.g., subsidence or uplift, east or west). Our results show a vertical mean velocity with a large percentage of negative values compared with positive values, with values ranging between −35.95 and 12.73 mm/year. On the other hand, horizontal mean velocity shows an equal percentage of negative and positive values, with values ranging from −22.98 to 21.85 mm/year. In the vertical mean velocity map, negative rates indicate subsidence and positive rates indicate land uplift, whereas in the horizontal mean velocity map, the negative values indicate deformations in the west direction and positive values indicate deformations in the east direction. It can be seen that the urban areas such as Cheria city are covered totally by negative values, which indicates that the whole city subsides (e.g., a sinkhole collapsed in 2009). However, rural areas with high-density borehole wells register the highest fluctuation of subsidence/uplift velocity rates (up to −35 mm/year) due to ongoing unsupervised excessive underground water pumping.
Five GCP points (P1, P2, P3, P4, and P5) were selected for the cumulative analysis of ground movement in the study area during the period between 2016 and 2022 (
Figure 6). The analysis of five ground control points illustrates a high rate of vertical movement caused by land subsidence. Overall, the PS-InSAR results are in close agreement with GPS measurements, with a ±3 mm/year difference. Point P1 shows an aggressive land subsidence at the start, reaching values greater than −275 mm, followed by a slow subsidence of approximately 30 mm registered at the end of the observation period. Point P2 (near the sinkhole that collapsed in 2004) shows that progressive land subsidence can reach a value of −200 mm. At P3, the progressive land subsidence reached −140 mm, and at P4, the uplift reached 25 mm, and huge land subsidence was registered at ~−250 mm. For P5, the progressive subsidence reached approximately −265 mm to −275 mm.
To analyze the subsidence trend in the study area, two profiles (AB and CD) with NE–SW and NWW–SEE directions were selected and plotted (see
Figure 7). The trend of cumulative vertical and horizontal ground movement profiles during the period of the analysis was visualized to determine the changes in the study area. The fluctuation of the vertical ground movement demonstrates that the study area suffers from a heavy subsidence pattern with a maximum displacement rate ranging between 0 and 350 mm, approximately. For the CD profile, the vertical ground movement shows a varying amount of subsidence where the maximum rates are between −100 and −500 mm.
An overview of eight sites selected from cumulative vertical ground movement profiles is presented in
Figure 8, with personalized zoom to monitor the total displacement resulting from over five years of analysis. For further detail, in
Figure 8b,c, the selected points A and B represent areas facing serious land subsidence highlighted in the map, with a sub-circular shape of the total value reaching more than −350 mm. In
Figure 8c,d, point C is situated near the urban perimeter and is represented on the map by near-zero vertical movements that indicate relatively stable zones. Point D suggests a medium to high land subsidence, with a value reaching −270 mm; while in
Figure 8e,f, points 1 and 4 refer to ground subsidence with values between −140 and −170 mm, respectively. On the other hand, points 2 and 3 represent areas with the highest land sinking values registered during the analysis period, at more than −490 mm. GPS distribution is very similar to the InSAR deformation field.
Our vertical velocities, obtained using the GPS rapid static mode, confirm the subsidence detected with InSAR.
The displacement rates derived from the PS-InSAR and GPS measurements were compared by calculating the root mean square error (RMSE) for vertical and horizontal displacement rates (
Table 2). The data are in close agreement, with an accuracy of 2–3 mm/year. RMSE comparison results show 2.8374 mm/year for vertical movement direction and 2.9155 mm/year for horizontal movement direction.
5. Discussion
Over the years, a limited focus on InSAR subsidence was given in Algeria, especially to subsidence-related issues in karstic terrains [
25]. Here, we compare the map of vertical ground deformation rates generated by PS-InSAR with GPS measurements. The results obtained in this research demonstrate the effectiveness of the monitoring system based on the integration of PS-InSAR data and GPS measurements. Our results (see
Figure 5) show that the study area is largely dominated by subsidence and needs continuous monitoring to understand the risks involved. We also noticed remarkable ground uplift at places with deformation ranging between 0 and +12.73 mm/yr. A possible explanation for this anomalous phenomenon could be related to the sudden hydrologic overpressure in the karst networks, seasonal water changes, and groundwater redistribution leading to the natural pressurization of the underground karst flow channels. These water storage losses correlate to flexural unloading [
26,
27]. Previous research has found that fluid extraction can lead to strain and an increase in pore pressure variations in the sediments [
28]. Slow land uplift has been proposed to be connected to increases in pore pressure accompanying groundwater redistribution [
29].
The reason behind the high subsidence rates detected in the study area seems to be a combination of multiple factors. The study area has known heavy groundwater exploitation in the last and present century. The main water resource for human supplies and agricultural uses in this area is groundwater. An in situ investigation of the consumption of underground water conducted in 2007 [
15] suggests that the average consumption is approximately 22 hm
3/year. The population grew by an average of 3.7%— from 74,129 in 2008 to 96,827 in 2020—in an area of about 267 km
2. Today, the study area suffers from groundwater drawdown due to two major factors: (a) the ongoing unsupervised groundwater exploitation to fulfill the increased water demand, and (b) the severe droughts and soil degradation affecting the area [
15]. In scientific literature (e.g., [
30,
31]), a world map climate classification was published which suggests that the study area is recognized as a hot desert climate type. Global warming in semi-arid regions makes drought more severe and last longer [
32], which makes dry regions drier [
15]. In the case of the study area, due to the lack of surface water resources needed for sustaining agriculture irrigation, local farmers rely on underground water as the only resource available. The study area is a semi-arid region that suffers from low rainfall precipitations and extreme heat variation. Therefore, there is huge reliance on underground water supplies for Cheria city and adjacent cities. Water pumped from underground mostly goes to: (a) agriculture which uses a huge amount of water in the form of primitive irrigation (see
Figure 1, below)—the water used is mostly lost to evaporation, and the remaining water either seeps underground to accelerate the dissolution of karst formation (see
Figure 2 and
Figure 3) and/or causes the swelling of clayey formations; and (b) the local drinking water network (Cheria city and adjacent cities), and subsequently to poorly maintained sewer systems (most Algerian cities do not have water treatment facilities) that usually dump water, such as irrigation water, on the surface.
Three boreholes (B1, B2, and B3; see
Figure 1 for location) were monitored between 2011 and 2012 (
Figure 9). The graphs in
Figure 9 show that the water table declined gradually without picking up. This comes with the knowledge that the current status of the underground water table is at its lowest and cannot fulfill human life needs for much longer. In the 1980s, the study area had a water table near the ground surface and even small swamps were documented [
15]. Excessive underground water pumping leads to water table decline, compression of soils, and land sinking. In semi-arid regions, the water table declines rack up rapidly. The drinking water crisis explodes with expanding water supply demand. The shortage in water halts agriculture-related activities. The consumption of huge quantities of the already dwindling water reserves creates an underground imbalance due to the extracted water densities that existed in the karst cavities. Therefore, soils covering these surfaces sink massively. This results in a subsidence phenomenon that will transform into sinkholes at any moment, especially in karstic formations. Furthermore, it was noticed that for a lot of sinkhole incidences in the region, the local authorities were not aware of them. Small subsurface circular areas of land sink gradually, and local farmers try to conceal them by adding soil on top in an attempt to stop land degradation and surface sinking. Additionally, local farmers and communities are afraid of land being seized without compensation. Algerian laws issued in 2004 concerning natural hazards do not recognize a lot of hazards such as climate change, heatwaves, land subsidence, etc. In addition to the lack of recognition, these laws were not updated to reflect the results of research reported by the scientific community. Therefore, no compensation was provided to landowners in the form of insurance or aid by local municipalities in the case of sinkhole collapse.
Based on the profiles of cumulative vertical ground movement fluctuation in AB and CD, urban areas are affected by continuous minimum rates of land deformation which alternate between slightly up and down variations. On the other hand, rural areas are extremely affected by maximum downward/upward movement deformation rates. This analysis clearly shows that subsidence dominated the study area during the investigation period, and the area surrounding the city experienced high to very high subsidence, while the interior of the city experienced low subsidence. In detail, the rural areas suffer due to dense illegal water boreholes being implemented.
Comparing the GPS measurements with the PS-InSAR results (
Table 2 and
Figure 6,
Figure 7 and
Figure 8) shows satisfying results. PS-InSAR deformation rates are in close agreement with GPS measurements covering the period from 2016 to 2022. The RMSE was less than 3 mm/year. Since RMSE can range from 0 to ∞, there is no absolute good or bad threshold. Therefore, it is recommended to assess RMSE based on the case study and the observation samples available. The obtained RMSE value is great considering the scale of analysis, the budget available, and time consumption. Studies (e.g., [
33,
34,
35,
36,
37,
38,
39,
40]) have obtained similar RMSEs to this study despite the difference in the processing approach. The obtained RMSEs between GPS and PS-InSAR deformation rates can be improved by enhancing the SAR processing algorithms and GPS measurements.
The comparison of GPS and PS-InSAR goes beyond comparing the deformation values from each method. Precise GPS positions can be used to enhance the deformation accuracy of InSAR images. PS-InSAR and GPS are two different technologies to monitor land surface deformation. Where PS-InSAR can measure changes in elevation over large areas with high accuracy, GPS can only measure elevation changes over small areas because it is time-consuming and requires high-cost equipment. Innovative methods of comparing and integrating GPS and InSAR measurements will facilitate enhanced land deformation mapping and provide a better understanding of subsidence processes.
PS-InSAR has been proven to be a very useful technique for monitoring deformations affecting rural and urban areas and provides information for large-scale and slow surface movement with centimeter to millimeter accuracy for hazards such as ground subsidence, structure collapse, landslides, mining subsidence, etc. Meanwhile, our results successfully demonstrated subsidence phenomena in the study area, but they can be further modified with in situ data analysis and other techniques, such as SBAS or Quasi-PS. It is also suggested that, in the future, a multi-scale (space-based and ground-based) study should be conducted to thoroughly analyze ground subsidence to avoid serious damage in this area. That being said, the applicability of the methodology presented here is subject to two conditions: (1) the availability of InSAR data covering the study region with high coherence; (2) the availability of measurements of some discrete variables highly correlated with InSAR data. In this study, we used GPS to assess the PS-InSAR results. Piezometric data can be also used, but the mathematical relationship with ground deformation is not established yet. The Sentinel-1 satellites, with a short revisit time, provide greater spatial coverage, and temporal sampling provides a great opportunity for developing countries to evaluate and monitor natural hazards with great consistency instead of GPS observation networks, which are very expensive and time-consuming. Our preliminary result from Sentinel-1 data demonstrates that excellent results were achieved over a wide area with the PS-InSAR, which is a key condition in order to apply the proposed methodology.