Quantifying Ground Subsidence Associated with Aquifer Overexploitation Using Space-Borne Radar Interferometry in Kabul, Afghanistan

Rapid population growth combined with recent drought events and decades of political instability have left the residents of Kabul facing water scarcity, significantly relying on groundwater. Groundwater overexploitation might have induced various magnitudes of ground subsidence, however, to date, no comprehensive study of ground subsidence in Kabul has been conducted. In this study, we investigated the spatio-temporal evolution of ground deformation phenomena and its main governing processes in Kabul from 2014 to 2019 using C-Band Sentinel-1 derived Interferometric Synthetic Aperture Radar (InSAR) time-series from both ascending and descending orbits to extract the two-dimensional (2D) surface displacement field. Four subsidence bowls were distinguished with highly variable spatial extents and deformation magnitudes over four separate aquifer basins, with the maximum value of −5.3 cm/year observed in the Upper Kabul aquifer basin. A wavelet analysis suggests that there is a strong correlation between the groundwater level variations and subsidence. Investigation of hydrogeological data further reveals that the observed subsidence could be attributed to the presence of highly compressible clayey soils. This detailed space-borne regional survey provides new insights into the main governing mechanism of land subsidence in Kabul and may direct better mitigation plans of potential hazards.


Introduction
Afghanistan is characterized by an arid climate with limited availability of surface water resources, hence adequate water supply is achieved only by drilling wells into aquifers to access groundwater [1,2]. Historically, the local population depended on water delivered from karez systems, which utilize groundwater resources from unconfined aquifers in alluvial fans replenished by snowmelt during spring months. Karezes, dating back 3000 years, are traditional well-like sub-horizontal tunnels that tap subterranean water to the surface by gravity [3,4]. As Afghanistan's population expands, there is an increasing need to withdraw more groundwater for domestic, agricultural, and municipal purposes. In particular, the capital city of Kabul, as one of the fastest growing and most water-stressed cities in the world, has been intensively exploiting aquifers, causing water levels to plummet. The main reason for this is a disturbed aquifer equilibrium with the groundwater discharge rate exceeding the recharge rate due to below-average rainfall periods and increased water consumption.
Aquifer overexploitation, altering the geodynamic state of the environment, might result in adverse environmental and socio-economic consequences such as reduction of a streamflow which, in turn, motion due to groundwater seasonal variations in the Po Plain, Italy was identified by Boni et al. [26] using the SqueeSAR algorithm applied to ERS-1/2 (from 1992 to 2000) and Radarsat (from 2003 to 2010) data. Other ground motion triggering factors in the area include clay shrink-swell process, load of new buildings, and tectonic uplift. Haghighi et al. [27] constructed a decadal time-series of deformation from multi-sensor observations in Tehran, Iran, with the results revealing the maximum subsidence rate exceeding 25 cm/year. The InSAR-derived time-series was compared with groundwater level data and was found to be the main cause of ground deformation in the area. Qu et al. [28] implemented the Stanford Method for Persistent Scatterers (StaMPS) to investigate land subsidence over Xi'an, China, using Envisat, ALOS, and TerraSAR-X during 2005-2012. They find a strong correlation between the ground subsidence profile and the depth to water table in the area. The mentioned studies benefit from time-series analysis that enables comparison of the InSAR-derived results with ancillary data (hydrogeological, geotechnical, etc.) in order to reveal either positive or negative relationships and recognize the driving mechanisms of the ground deformation.
To date, no comprehensive study of land subsidence in Kabul has been conducted. In this study, for the first time, we mapped and quantified the ground surface deformation over the study area using C-band Sentinel-1 SAR products during the period from 2014 to 2019. Then, we implemented spatio-temporal analysis of ground motion to infer the triggering factors and demonstrate that InSAR-derived ground deformation rates have patterns consistent with observed groundwater changes. We also established the direct relationship of geotechnical properties of the area with the magnitude of subsidence. The results have considerable implications for developing more sustainable water management strategies and efficient groundwater pumping or artificial replenishing schemes for the region.
The paper is structured as follows: Section 2 describes the geographical and geological settings of the study area and presents the problem statement; SAR data and InSAR processing chain employed are outlined in Section 3; Section 4 presents the InSAR-derived results, whilst Section 5 discusses and analyses the findings and factors triggering ground subsidence. The main conclusions are summarized in Section 6.

Descriptions of the Study Area
Kabul is the capital and largest city in Afghanistan with a population of 4.2 million. It is located in the southernmost part of the Kabul basin, bounded by Hindu Kush mountain ranges-Paghman mountains to the west as high as 4400 m, and the Kohe Safi mountains to the east reaching 3000 m in altitude. The city is divided into upper and lower parts by the Guzarga-Asmai mountain. The mountains are made up of the Precambrian metamorphic rocks, including gneiss (Xgn) and schist (Ym) (Figure 1b). The southern rims of the city are dominated by mountains of younger age (upper Paleozoic-Mezozoic) such as limestone, dolomite, sandstone, and siltstone (Cambrian or Neoproterozoic limestone and dolomite (ZeEld), late or Middle Triassic sandstone and siltstone (T23ssl), middle or early Triassic limestone, dolomite and marl (T12ld), silurian limestone and marl (Sls), late Permian limestoneand dolomite (P2ld)) [29]. Lowlands and valleys of the basin are deposited with unconsolidated Tertiary (Neogene) and Quaternary alluvial, colluvial, and lacustrine sediments as a result of erosion. Sediments are generally coarse-grained near the mountain foothills, such as gravel and pebble-sized deposits from alluvial fans (Q34ac), with fining trend to the center of the valleys leading to fine-grained sand and silt. The youngest formation of the basin is loess (Q3loe) of Pleistocene epoch [30,31]. Altitudes of the central plains range from 1800 m in Lower Kabul and Logar, and 2200 m in Paghman and Upper Kabul [32].
Three major rivers flow through the city: Kabul, Logar, and Paghman. Upper Kabul is drained by upper Kabul and Paghman rivers, that join to form lower Kabul river as it passes through a narrow gorge in the Guzarga-Asmai mountain. The lower Kabul is drained by Kabul and Logar rivers that also join to form a single channel-Kabul river and pass through yet another mountain gorge [33]. The city is underlain by four interconnected aquifers-Upper Kabul, Lower Kabul, Logar, and Paghman-each Blue and red rectangles outline the footprints of ascending and descending Sentinel-1 SAR datasets indicated as TA71 and TD78, respectively. The grey frame shows the area of reconstructed 2D surface deformation field. (b) The geological map of the study area corresponding to the grey frame [36].
Since 1960, drought and overexploitation of groundwater for multiple competing purposes led to progressive groundwater table drawdown and depleted aquifers in Kabul basin. As a result, many of the shallow wells, springs, and karezes have dried up. Banks et al. [30] reported declines in the water table elevation of 4 to 6 m in the Kabul area as a result of the drought in 1998-2002. Although groundwater levels in some areas of the Kabul basin outside the urban area have been rising as a response to a rainfall increase back to the level prior to the drought since 2004, the water table within the city boundaries has continued to decrease steadily due to effects of overabstraction superimposed upon climatic trends [37]. Groundwater data obtained from the National Water Affairs Regulation Authority (NWARA) of Afghanistan shows the borehole yields were dramatically reduced by up tõ 20 m in 5 years from 2014 to 2019 in Upper and Lower Kabul sub-basins, as illustrated in Figure 2. It is evident from the plot that the largest drop occurred in Upper Kabul basin, while the water level in Logar shows seasonal variations without significant decreasing trend. ogressive groundwater table drawdown and depleted aquifers in Kabul basin. As a r of the shallow wells, springs, and karezes have dried up. Banks et al. [30] reported declin ater table elevation of 4 to 6 m in the Kabul area as a result of the drought in 1998ugh groundwater levels in some areas of the Kabul basin outside the urban area have as a response to a rainfall increase back to the level prior to the drought since 2004, the w within the city boundaries has continued to decrease steadily due to effects of overabstra imposed upon climatic trends [37]. Groundwater data obtained from the National W rs Regulation Authority (NWARA) of Afghanistan shows the borehole yields atically reduced by up to ~20 m in 5 years from 2014 to 2019 in Upper and Lower K asins, as illustrated in Figure 2. It is evident from the plot that the largest drop occurr r Kabul basin, while the water level in Logar shows seasonal variations without signif asing trend.  Figure 5. Precipitation time-series is retrieved from CRU TS 4.03 [38].

ta and Methods
Sentinel-1, a space mission of the European Space Agency (ESA), has resulted in a para in the field of InSAR: configurations of the constellation, small perpendicular and tem ines, and free data availability [39] allow accumulation of a long-term archive of glo stent data. A stack of 106 ascending and 101 descending Sentinel-1 Single Look Complex ( ucts obtained in Interferometric Wide Swath (IWS) mode spanning the period from Oc to May 2019 were utilized in this study. SAR data parameters used for this study are desc ble 1, while the spatial and temporal baselines of the generated interferograms are show e 3.

Data and Methods
Sentinel-1, a space mission of the European Space Agency (ESA), has resulted in a paradigm shift in the field of InSAR: configurations of the constellation, small perpendicular and temporal baselines, and free data availability [39] allow accumulation of a long-term archive of globally consistent data. A stack of 106 ascending and 101 descending Sentinel-1 Single Look Complex (SLC) products obtained in Interferometric Wide Swath (IWS) mode spanning the period from October 2014 to May 2019 were utilized in this study. SAR data parameters used for this study are described in Table 1, while the spatial and temporal baselines of the generated interferograms are shown in Figure 3.

Interferometric Processing
Interferometric processing was implemented using the GAMMA software [40]. Following the data import, radiometric calibration, and mosaicking of the bursts and sub-swaths, all of the SLC SAR images were co-registered and resampled to a reference SAR image to produce the stack of interferograms using the precise orbit information from ESA. Interferometry implies pixel-to-pixel and finer subpixel-to-subpixel matches between common features in image pairs. Terrain Observation by Progressive Scans (TOPS) interferometry particularly requires highly precise alignment and co-registration in the azimuth direction because of the signal properties resulting from the antenna beam steering from aft to the fore [41]. Phase jumps between adjacent bursts resulted from the strong along-track Doppler Centroid variations are observed unless co-registration accuracy of a few thousands of a pixel is achieved [42]. To reach a high level of co-registration accuracy, resampling using an intensity matching procedure followed by refinement using spectral diversity approach that considers the interferometric phase of the subsequent bursts' overlap area were applied iteratively [43]. An azimuth spectrum deramping technique is applied to eliminate along-track Doppler Centroid which runs through a spectral ramp.
An interferogram contains information including displacement (dynamic), topography (static), and errors due to atmospheric disturbance, orbit, and decorrelation noise. All of the phase components other than that due to displacement should be removed from the interferograms using temporal-spatial filtering. In this study, the external Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) with 3-arcsecond (90 m) resolution was used to remove the topographic component and geocode the differential interferograms. The flat-earth phase was removed using the precise orbit data of the satellites from ESA. To enhance the signal-to-noise ratio of the interferograms and aid the following unwrapping procedure, a multi-looking operation with two looks in the azimuth direction and ten looks in the range direction was employed. Phase information of differential interferograms was sampled in a discrete wrapped phase in the interval (−π, +π]. The phase unwrapping procedure must be applied to relate the interferometric phase of the wave to absolute distance, thus an integer of 2π must be added to derive the unambiguous phase

Interferometric Processing
Interferometric processing was implemented using the GAMMA software [40]. Following the data import, radiometric calibration, and mosaicking of the bursts and sub-swaths, all of the SLC SAR images were co-registered and resampled to a reference SAR image to produce the stack of interferograms using the precise orbit information from ESA. Interferometry implies pixel-to-pixel and finer subpixel-to-subpixel matches between common features in image pairs. Terrain Observation by Progressive Scans (TOPS) interferometry particularly requires highly precise alignment and co-registration in the azimuth direction because of the signal properties resulting from the antenna beam steering from aft to the fore [41]. Phase jumps between adjacent bursts resulted from the strong along-track Doppler Centroid variations are observed unless co-registration accuracy of a few thousands of a pixel is achieved [42]. To reach a high level of co-registration accuracy, resampling using an intensity matching procedure followed by refinement using spectral diversity approach that considers the interferometric phase of the subsequent bursts' overlap area were applied iteratively [43]. An azimuth spectrum deramping technique is applied to eliminate along-track Doppler Centroid which runs through a spectral ramp.
An interferogram contains information including displacement (dynamic), topography (static), and errors due to atmospheric disturbance, orbit, and decorrelation noise. All of the phase components other than that due to displacement should be removed from the interferograms using temporal-spatial filtering. In this study, the external Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) with 3-arcsecond (90 m) resolution was used to remove the topographic component and geocode the differential interferograms. The flat-earth phase was removed using the precise orbit data of the satellites from ESA. To enhance the signal-to-noise ratio of the interferograms and aid the following unwrapping procedure, a multi-looking operation with two looks in the azimuth direction and ten looks in the range direction was employed. Phase information of differential interferograms was sampled in a discrete wrapped phase in the interval (−π, +π]. The phase unwrapping procedure must be applied to relate the interferometric phase of the wave to absolute distance, thus an integer of 2π must be added to derive the unambiguous phase difference [44]. The minimum cost flow method with the coherence threshold set to 0.4 was adopted to unwrap the differential interferograms [45,46].

Time-Series Analysis
Following the unwrapping of the differential interferograms, the time-series analysis was performed. A variety of time-series algorithms have been developed in recent decades, such as Persistent Scatterers, Small Baseline Subset, and various hybrid methods such as SqueeSAR [47][48][49]. The SBAS technique uses a number of multi-prime differential interferograms with small orbital separation and temporal intervals to minimize the spatial decorrelation phenomena. In this study, we used a modified SBAS method, which integrates the atmospheric artifact estimation model, referred to as InSAR time-series + atmospheric estimation model (InSAR TS + AEM). The phase contribution resulting from temporal atmospheric delay between subsequent observations must be calculated and removed. Taking into account the fact that the atmospheric signal is highly correlated in space but not in time [50], whereas non-linear displacements are correlated both in space and time, a spatial low-pass filter together with a temporal high-pass filter was used to distinguish atmospheric effects from deformation signal. The Atmospheric Phase Screen (APS) is estimated using a temporarily linear velocity (TLV) model and filtered out from the unwrapped phase. Time-series analysis utilizes partially coherent pixels together with fully coherent pixels over time. Displacements measured by InSAR are given relative to a reference point and a reference acquisition date. The reference point was selected to the north-west of the city, based on the following criteria: (1) situated away from the urban area and (2) exhibits high coherence values (longitude and latitude at 69.0 and 34.6, respectively). The interferogram network was inverted with respect to the first SAR acquisition date. The time-series approach implemented in the study is based on the methodology proposed by Li et al. [13].

Multi-Geometry Data Fusion to Derive 2D Model
The traditional InSAR technique is capable of measuring Line-Of-Sight (LOS) displacement only-the vector in the satellite-to-ground direction-whereas real deformation vectors occur in vertical (up) and horizontal (east-west, north-south) directions. In areas with overlapping spatiotemporal coverage, multi-geometry InSAR data fusion can be utilized to discriminate 2D or 3D components of displacement. Therefore, in order to overcome the limitation of one-dimensional LOS measurements and yield reliable solution, we combine two independent LOS displacements observed along the LOS from the satellite in ascending and descending orbits to retrieve vertical and horizontal components of deformation. At least three independent LOS measurements are required to solve for the 3D ground deformation pattern [51]. There are no independent datasets such as GNSS available to integrate with LOS measurements in order to derive all three components in Kabul. Considering near polar orbits and the incidence angle of current SAR satellites, InSAR observations are mostly sensitive to the vertical component while their sensitivity to the north component is minimal, therefore it is assumed that the detected displacement is caused by vertical and E-W motion, and the N-S motion is neglected [52]. Other techniques such as combination of multi-pass LOS and azimuth measurements from Pixel Offset Tracking or Multiple Aperture Interferometry (MAI) were also not implemented because the magnitude of the land subsidence events is not large enough to derive adequate results-the techniques mentioned above are valid for large motion events at a meter scale [21,53].
The equation system Equation (1) is solved to obtain the vertical and E-W components of the displacement from the ascending and descending LOS motions [54]: where u los is the LOS displacement detected for each geometry, θ asc and θ dsc are the incidence angle determined for each pixel, and α asc and α dsc are the heading angles of the satellite, positive clockwise from the north.

Wavelet Transform Analysis
A signal decomposition algorithm, widely used for time-series analysis, has been implemented in order to acquire insights about the physical relationships between groundwater level and ground deformation. It enables determining whether the two datasets are suggestive of causality by decomposing a time-series over a time-frequency space, identifying a common power, and characterizing their phase relationship. Wavelet transform expands on the classical Fourier transform, but allows better temporal resolution [55]. First, Continuous Wavelet Transform (CWT) examines the evolution of each phenomenon over time by detecting the dominant oscillating modes. To do so, a wavelet is applied as a bandpass filter to the input signal. This study employed the Morlet wavelet, which consists of a sine wave modulated by a Gaussian window. CWT is followed by the cross wavelet transform (XWT), which measures the covariance of two datasets by convolution of the CWT-derived scalograms of two signals resulting in a cross-scalogram. XWT detects regions with high common time-localized oscillations. Wavelet coherence (WTC) estimates the coherence of the two series even though the common power is low, thus revealing the locally phase-locked correlation. A cause-effect relationship is inferred for patches with common high powers on XWT and WTC with the use of phase arrows revealing either in-or anti-phase relationship. The significance level is determined using Monte Carlo methods [56,57]. Cross-correlation of datasets allows computing the time lag between the cause-effect events using the following equation (Equation (2)): where ∆ ϕ is the phase difference and T is the period [58]. Data pre-processing includes the following steps: (1) ensuring that irregularly-spaced datasets have a uniform time increment by interpolating missing gaps between the consecutive InSAR acquisition dates; (2) interpolating groundwater level time-series to a more temporally dense InSAR time-series; and (3) removing the trend from InSAR-derived results. Four subsiding bowls with spatially heterogeneous pattern and temporally varying rates were delineated in LOS displacement maps of both geometries, accounting for the maximum mean annual rates of −4.08 and −4.3 cm/year for ascending and descending tracks, respectively, during the study period from October 2014 to May 2019. The maximum cumulative displacements account for −17.9 and −19.6 cm for ascending and descending tracks, respectively. The deformation fields retrieved from both geometries show similar spatial patterns. To quantitatively validate the agreement between two tracks' results, correlation analysis was implemented, as illustrated in Figure 4e. The high correlation coefficient of 0.75 and the low Root Mean Square Error (RMSE) of 0.60 cm/year suggest that the two measurements are consistent. The good agreement of two independent LOS displacement values implies that vertical deformation is dominant in the area with limited atmospheric phase contribution after the APS filter. Observed minor differences should be attributed to different looking angles of geometries.  Four subsiding bowls with spatially heterogeneous pattern and temporally varying rates were delineated in LOS displacement maps of both geometries, accounting for the maximum mean annual rates of −4.08 and −4.3 cm/year for ascending and descending tracks, respectively, during the study period from October 2014 to May 2019. The maximum cumulative displacements account for −17.9 and −19.6 cm for ascending and descending tracks, respectively. The deformation fields retrieved from both geometries show similar spatial patterns. To quantitatively validate the agreement between two tracks' results, correlation analysis was implemented, as illustrated in Figure 4e. The high correlation coefficient of 0.75 and the low Root Mean Square Error (RMSE) of 0.60 cm/year suggest that the two measurements are consistent. The good agreement of two independent LOS displacement values implies that vertical deformation is dominant in the area with limited atmospheric phase contribution after the APS filter. Observed minor differences should be attributed to different looking angles of geometries.

Decomposition into Horizontal and Vertical Displacements
The spatial distribution of the subsidence bowls obtained after decomposition into horizontal and vertical components is shown in Figure 4c

Decomposition into Horizontal and Vertical Displacements
The spatial distribution of the subsidence bowls obtained after decomposition into horizontal and vertical components is shown in Figure 4c

Spatial Patterns
Four subsidence bowls with spatially non-uniform patterns were detected in all four aquifers: Upper Kabul Basin (UKB), Paghman, Lower Kabul Basin (LKB), and Logar, with the maximum cumulative subsidence at each reaching 23.8, 10.9, 12.9, and 8.6 cm, respectively (Figure 5a,b). The Guzarga-Asmai mountain, extending from south to north in the middle of the city and splitting it into two, acts as a barrier between the deforming areas. In unconsolidated sedimentary deposits, an increase in effective stress results in compaction as opposed to incompressible consolidated metamorphic rocks. In addition, areas along the rivers flowing through the center of the city are also less exposed to deformation phenomena. Particularly looking closer at areas along the Paghman river on the subsidence map, one can note a distinct non-deforming narrow strip of ground between two highly deforming bowls. This may be explained by relative "saturation" of the soil proximal to the river as a result of lateral movement of water from areas of greater hydrostatic pressure to those of lesser [60], and therefore increased pore water pressure and decreased effective stress. On the other hand, the surface geology of the river banks (Figure 5c) comprised of coarse conglomerate and sandstone (Q4a) may also contribute to their relative stability in contrast with adjacent layers of soft compressible loess (Q3loe). Thus, a mountain extending from south to north, and rivers running from east to west, largely determine the spatial extent of deformation, confining the subsidence zones and resulting in four spatially discontinuous subsidence bowls.
In order to quantitatively assess the ratio of subsidence occurrence in each of the above-mentioned units-river banks, mountain, and unconsolidated sedimentary The Guzarga-Asmai mountain, extending from south to north in the middle of the city and splitting it into two, acts as a barrier between the deforming areas. In unconsolidated sedimentary deposits, an increase in effective stress results in compaction as opposed to incompressible consolidated metamorphic rocks. In addition, areas along the rivers flowing through the center of the city are also less exposed to deformation phenomena. Particularly looking closer at areas along the Paghman river on the subsidence map, one can note a distinct non-deforming narrow strip of ground between two highly deforming bowls. This may be explained by relative "saturation" of the soil proximal to the river as a result of lateral movement of water from areas of greater hydrostatic pressure to those of lesser [60], and therefore increased pore water pressure and decreased effective stress. On the other hand, the surface geology of the river banks (Figure 5c) comprised of coarse conglomerate and sandstone (Q4a) may also contribute to their relative stability in contrast with adjacent layers of soft compressible loess (Q3loe). Thus, a mountain extending from south to north, and rivers running from east to west, largely determine the spatial extent of deformation, confining the subsidence zones and resulting in four spatially discontinuous subsidence bowls.
In order to quantitatively assess the ratio of subsidence occurrence in each of the above-mentioned units-river banks, mountain, and unconsolidated sedimentary deposits-Frequency Ratio analysis was performed. The Frequency Ratio (FR) model is a bivariate statistical model designed to evaluate the correlation between the conditioning factor and the phenomenon. It was calculated following the methodology outlined by Pradhan et al. and Bianchini et al. [61,62]. A FR ratio higher than 1 means high a correlation with subsidence, while FR value lower than 1 indicates a lower correlation. According to the results, the FR for the river banks, mountain, and unconsolidated sediments account for 0.53, 0.0, and 3.9, respectively, and confirm the positive relationship between the geological formations and subsidence distribution.
Another pattern arising from the results is the presence of the zone of low deformation in Logar basin near well CW2 (Figure 5d). Two factors were found to be responsible for the phenomena: undisturbed groundwater flow evident from records of well CW2 ( Figure 6) and the lithological setting of the area. As illustrated in Figure 5d, surface geologic units along the river include Quaternary and Neogene coarse conglomerate and sandstone (Q2a and N2cgs), which is less compressible compared to loess (Q3loe).
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 21 deposits-Frequency Ratio analysis was performed. The Frequency Ratio (FR) model is a bivariate statistical model designed to evaluate the correlation between the conditioning factor and the phenomenon. It was calculated following the methodology outlined by Pradhan et al. and Bianchini et al. [61,62]. A FR ratio higher than 1 means high a correlation with subsidence, while FR value lower than 1 indicates a lower correlation. According to the results, the FR for the river banks, mountain, and unconsolidated sediments account for 0.53, 0.0, and 3.9, respectively, and confirm the positive relationship between the geological formations and subsidence distribution. Another pattern arising from the results is the presence of the zone of low deformation in Logar basin near well CW2 (Figure 5d). Two factors were found to be responsible for the phenomena: undisturbed groundwater flow evident from records of well CW2 ( Figure 6) and the lithological setting of the area. As illustrated in Figure 5d, surface geologic units along the river include Quaternary and Neogene coarse conglomerate and sandstone (Q2a and N2cgs), which is less compressible compared to loess (Q3loe). To examine whether the subsidence spatial distribution is controlled by local geological faulting, Quaternary surface fault traces obtained from [59] were superimposed on the cumulative subsidence map, as illustrated in Figure 5a. Most of the geological structures appear to be absent in the city boundaries and cannot be held to be responsible for observed deformation phenomena in Kabul. However, it is still possible that some of the deformation may be related to the unmapped minor faulting system around the city. In addition, deformation around the UKB may have contributed considerable strain accumulation of faults within Kabul block and increased the chance of future seismic events.

Temporal Patterns
Although the annual mean velocity and cumulative values may be valuable to recognize processes with linear trends, they may not be effective to detect non-linear motion, accelerationdeceleration events. and rapid displacements depending on the employed unwrapping approach [58]. To investigate the temporal evolution of highly subsiding areas, the time-series analysis of To examine whether the subsidence spatial distribution is controlled by local geological faulting, Quaternary surface fault traces obtained from [59] were superimposed on the cumulative subsidence map, as illustrated in Figure 5a. Most of the geological structures appear to be absent in the city boundaries and cannot be held to be responsible for observed deformation phenomena in Kabul. However, it is still possible that some of the deformation may be related to the unmapped minor faulting system around the city. In addition, deformation around the UKB may have contributed considerable strain accumulation of faults within Kabul block and increased the chance of future seismic events.

Temporal Patterns
Although the annual mean velocity and cumulative values may be valuable to recognize processes with linear trends, they may not be effective to detect non-linear motion, acceleration-deceleration events. and rapid displacements depending on the employed unwrapping approach [58]. To investigate the temporal evolution of highly subsiding areas, the time-series analysis of pixels was adopted. The location of pixels was selected from all four aquifers to represent the highest magnitude in each. As illustrated in Figure 5b, the time-series of vertical deformation generally exhibit progressive evolution pattern for all deforming areas yet at different rates. While highly subsiding UKB exhibits a linear and steady trend, the pattern in Logar oscillates to a small degree, which may indicate the elastic behavior of deformation in the area. The rate of subsidence at UKB is the largest in the city, twice and three times larger than the other sites.

Discussion
A wide range of factors govern the extent and intensity of ground deformation, including groundwater pumping, compaction of unconsolidated sediments, rapid urbanization, and tectonic activity. Therefore, the InSAR-derived subsidence results were interpreted in relation to the spatio-temporal variations of hydrological and geotechnical conditions of the study area to identify the governing factors of the deformation. In this section, we introduce the insights emerging from the results and potential mechanisms inducing the subsidence.

Comparisons with Groundwater Level Change
For any arbitrary surface below the groundwater table, pore water pressure (PWP) and effective stress balances the total stress from the overlying strata [63]. Excessive pumping of groundwater leads to decreased PWP and increased effective stress in unconsolidated sediments that puts overburden on the soil matrix from overlying weight and causes consolidation of the soil, which is manifested at ground surface as subsidence. This is explained by the one-dimensional consolidation theory by Terzaghi, which states that compression of soils occurs when the soil voids are dewatered and the weight is transferred from pore water to the granular structure of soil [7,64].
There is a progressive groundwater level drawdown in the majority of wells in UKB, LKB, and Paghman aquifers, whereas the water level is predominantly stable in Logar, showing some seasonal fluctuations: the peak is observed in April-May and the minimum in October-November. The largest drawdown is observed in the Upper Kabul aquifer. To investigate the relationship between temporal variations of ground subsidence and groundwater level changes in Kabul, displacement time-series were compared against the monthly standing water level (SWL) ( Figure 6). SWL data from 2014 to 2019 were analyzed to match the InSAR observation time interval. Well locations were selected to represent differing deformation rates (low, moderate, high) in each basin, with three taken from each.
There is a consistent pattern emerging from the comparisons of deformation time-series and SWL: the subsidence pattern corresponds with the downward trend in the groundwater level, while small ground uplift-subsidence movements were detected in areas with seasonal fluctuations but no long-term groundwater drawdown. SWL at wells 25, 23, and 88 show nearly sinusoidal seasonal oscillations, as does the ground subsidence, showing patterns that fluctuate by a fraction of cm, as opposed to steady and linear subsidence patterns in areas experiencing long-term subsidence (e.g., well 80). The two contrasting cases are juxtaposed in the last two subplots of Figure 6. In addition, the magnitude of the ground surface deformation and volume of SWL change appear to reflect each other, for example, well 80; as the groundwater level decreased significantly by 29 m, the ground surface subsided considerably by 8.2 cm. Conversely, as the groundwater level dropped to a relatively small degree by 3 m at well 28, the cumulative ground subsidence accounted for 3.3 cm. These findings suggest that subsidence is largely a function of a groundwater level change in Kabul.
To further elaborate how the two time-series evolved against each other and reveal the relative lag time, wavelet transform analysis was implemented. Data pre-processing included the following steps: (1) ensuring that irregularly-spaced datasets have uniform time increment by interpolating missing gaps between the consecutive acquisition dates; (2) interpolating groundwater level time-series to a more temporally dense InSAR time-series; (3) removing trend from InSAR-derived results. The majority of subsidence in Kabul exhibited a linear trend, particularly in UKB, LKB, and Paghman basins due to the heavy pumping of groundwater. However, wavelet analysis is designed to deal with localized oscillations. Unlike the three basins mentioned above, the ground movement in Logar exhibits some fluctuations enabling implementation of wavelet transform analysis. The results of CWT, XWT, and WTC analysis at four wells-CW2, 25, 23 and 88-are shown in Figure 7a- Regarding well 23 (Figure 7c), there is no periodic signal presenting high power above 5% significance level in XWT, however WTC reveals two positively correlated patches at bands 4 and 8, corresponding to 74 and 148 days, respectively, from 2016 to 2018. In particular, around 2018 a region within a 95% confidence interval and coherence value of 0.8, in which direction arrows point right and up, evidences that ground subsidence follows a groundwater level drop and vice versa, The continuous wavelet transform shows the distribution of a frequency component with respect to time. The thick contour indicates the 5% significance level against red noise, while the area demarcated with the thinner black line is the cone of influence (COI). The COI suffers from edge effects of signal analysis therefore does not yield reliable results. Warm colors denote high power whereas low power is indicated by cold colors. XWT and WTC can reveal common high power of two signals in terms of correlation and coherence, respectively. The direction of arrows indicates the relative phase relationship: pointing right means in-phase (positively correlated); pointing left means anti-phase (negatively correlated); straight down implies that ground deformation (GD) leads groundwater standing water level (SWL) by 90 • ; and straight up, vice versa, implies GD lags behind SWL by 90 • .
According to the CWT power spectra of GD, due to the very low magnitudes of ground motion oscillations in Paghman basin, there are no high wavelet powers revealed at the four sites. SWL scalograms, however, reveal significant periodicities with power 8, which manifest throughout the whole study period from 2014 to 2019 at around band 20, which corresponds to a one-year period (370 days) and confirms seasonal groundwater level changes at all four sites.
Results of wavelet analysis for well CW2 are illustrated in Figure 7a. Seasonality of SWL is occasionally disturbed by groundwater pumping, hence scalogram exhibits high-power around band 20 (370 days) all the way to band 4 (74 days), reflecting groundwater level changes due to pumping at lower frequencies. XWT reveals a patch within a 95% confidence level yet with a low power between 2016 and 2017, direction arrows showing a shift by 310 • between two time-series. Nevertheless, WTC reveals regions with significant coherence despite the common power being low, which should be considered as a "local correlation" between the two time-series in the time-frequency space [56]. Unlike XWT, there is no requirement for large values of individual CWTs to produce high coherence; signal similarities in WTC depend only on the time-frequency pattern and not on the power of components in the CWT [65]. Three regions with high power are evident in WTC at bands 3, 4 and 8 with the phase arrows pointing straight down, up at 45 • , and rightward, respectively. The former implies SWL lagged behind GD at a frequency of 55 days in 2017 by 27 days, then SWL led GD every 74 days with the time lag of 9 days, and the latter implies an in-phase relationship.
According to the results of well 25 (Figure 7b), the XWT spectrum reveals a low magnitude region between 2018 and 2019 yet within the 5% significance level, with the phase difference of 230 • . The WTC reveals two discrete high-magnitude regions at bands 4 (74 days) and 6 (111 days). The former occurred in Autumn 2015, revealing am in-phase relationship between SWL and GD, while the latter in 2018 exhibits an anti-phase relationship. There might be two reasons for the negative correlation in 2018-either the presence of another factor that triggered ground deformation in 2018 or poor quality of raw SWL data.
Regarding well 23 (Figure 7c), there is no periodic signal presenting high power above 5% significance level in XWT, however WTC reveals two positively correlated patches at bands 4 and 8, corresponding to 74 and 148 days, respectively, from 2016 to 2018. In particular, around 2018 a region within a 95% confidence interval and coherence value of 0.8, in which direction arrows point right and up, evidences that ground subsidence follows a groundwater level drop and vice versa, groundwater level increase leads to uplift, at a frequency of~5 month. This suggests elastic behavior of deformation. However, taking into account the short time delay (16 days) between the events, we speculate that these high-frequency low-magnitude subsidence-rebound episodes at well 23 might rather result from seasonal shallow clay shrink-swell processes superimposed with long-term subsidence.
Results of wavelet analysis for well 88 are shown in Figure 7d. CWT of SWL showing an almost uniform shape of the high-power patch (band 20) confirms the largely undisturbed groundwater condition, where seasonal recharge and discharge are in equilibrium. Although XWT and WTC results are below 5% significance level against red noise. WTC shows relatively high-energy regions with phase arrows indicating 45 • between the two datasets at bands 3 and 6 between 2018 and 2019, revealing time lags of 7 and 13 days, respectively. We again interpret these high-frequency low-magnitude signals as clay shrink-swell episodes.
Results of wavelet transform analysis confirm the capabilities of the method to recognize linkages between two time-series and reveal their time lag. However, as ground subsidence is a complex phenomenon, often involving multiple driving factors, it is necessary to incorporate other auxiliary datasets to infer the relationships that wavelet analysis fails to capture or explain, such as with the case of well 25 in our study.

Comparisons with Geotechnical Properties of Sediments
The compressibility of soil is largely caused by the prevalence of fine-grained sediments in subsurface, particularly clay. The amount of water held in soil is linearly related to their clay content [66]; clay-rich sediments lose more water under decreasing moisture conditions, therefore to a greater extent are susceptible to shrinking compared to other soil types [67]. Given the growing urbanization of Kabul and associated increase in static external load from heavy buildings, compressible clay deposits are prone to compaction and subsidence even more. Conversely, the presence of sand particles in sediments reduces the consolidation rate of the soil because they reduce the amount of water held by the soil [68]. Assuming the same groundwater withdrawal rates across the network of wells, differential consolidation is a function of variation in thickness and type of lithology.
Taking into consideration the compaction properties of soils, ground deformation spatial patterns were compared against the lateral variation of superficial geology using available stratigraphy logs. To do so, subsidence profiles along two cross-sections in UPK and LKB were extracted (Figure 8). Soil stratigraphy is plotted against the SWL as depth to water level from the ground surface. It should be noted that sudden significant water level drops at well TW4 (autumn 2016 and 2017) correspond to pumped water level (PWL).
The examination of the cross-section B-B' reveals a strong impact of the lateral variation of the sediment's lithology on subsidence. Subsidence at all three wells show homogenous rates accounting for cumulative~4 cm subsidence, yet SWL decrease at TW1 and TW2 is twice as much as at TW3. This is due to the fact that shallow subsurface at well TW3 comprises 7 m of clay, as opposed to a thick layer of gravel and sand with some interbedded clay at the two other sites TW1 and TW2. Being at the two opposite extremes in terms of their compressibility characteristics, clays are considered highly compressible, while coarse-grained gravel and sand are the least. Thus, spatial variability of clay largely determines the spatio-temporal characteristics of subsidence in LKB. Shallow subsurface at cross-section A-A' at both wells exhibits similar subsurface asset comprising 6-7 m of clay, underlain with coarse-grained deposits such as gravel, sand, and pebble. Nonetheless, subsidence rates differ considerably at the two sites, accounting for 3.5 and 13 cm at wells TW4 and 6, respectively. Low subsidence rates at TW4 are explained by the relatively recent drop of groundwater level and drainage of water from sediment particles. The opposite is true about well 6, which experienced a significant long-term SWL decrease. The water table difference between the two accounts for 28 m. Due to the low permeability of clays, the dissipation of the excess pore pressure to reach the equilibrium will appear slowly, compared to other soil types. Accordingly, consolidation and subsidence at well TW4 is expected to take place following some time lag. Moreover, the consolidation rate of clays is time-dependent and it is important to determine whether primary or secondary consolidation is involved. Primary consolidation, occurring while the excess PWP is being dissipated, induces greater magnitudes of subsidence in a shorter period, while secondary consolidation, taking place due to plastic deformation and readjustment of soil particles, appears slowly over a longer period of time [69]. If the subsidence occurring at well TW4 involves the secondary consolidation, it will manifest even slower.
Moreover, the consolidation rate of clays is time-dependent and it is important to determine whether primary or secondary consolidation is involved. Primary consolidation, occurring while the excess PWP is being dissipated, induces greater magnitudes of subsidence in a shorter period, while secondary consolidation, taking place due to plastic deformation and readjustment of soil particles, appears slowly over a longer period of time [69]. If the subsidence occurring at well TW4 involves the secondary consolidation, it will manifest even slower. The presence of fine-grained sediments such as clay and silt is found to be the main culprit behind ground subsidence globally, which is accelerated by groundwater overdraft [70]. Our results are consistent with the findings of previous studies quantifying subsidence of megacities sitting on thick layers of compressible deposits due to reduced PWP, notable examples of which include Shanghai, Mexico, and Bangkok [71][72][73]. The recent rapid urbanization of Kabul is also likely to have The presence of fine-grained sediments such as clay and silt is found to be the main culprit behind ground subsidence globally, which is accelerated by groundwater overdraft [70]. Our results are consistent with the findings of previous studies quantifying subsidence of megacities sitting on thick layers of compressible deposits due to reduced PWP, notable examples of which include Shanghai, Mexico, and Bangkok [71][72][73]. The recent rapid urbanization of Kabul is also likely to have contributed to the ground subsidence, as the loading of buildings is transferred to the soil matrix. In addition to its high compressibility, the platy grain shape makes subsidence in clays inelastic (irreversible). Upon the original deposition of clay, grains are arranged in erratic orientations, thus having a larger volume and enough room to hold water in pores. However, once clays are drained, grains are rearranged as stacks, with a loose volume and are permanently compact [70]. This process of soil particle readjustment is part of the secondary consolidation, which might take up to a decade and last even after the groundwater withdrawal rates have been made sustainable.

Conclusions
In this study, multi-geometry Sentinel-1 InSAR data stacks acquired during 2014-2019 were used to infer the distribution, extent, and magnitude of the ground subsidence in Kabul due to groundwater overexploitation. Multi-geometry data fusion was implemented to extract vertical and horizontal components of the motion. The conclusions of the study are summarized as follows: 1.
The integration of ascending and descending LOS displacement rates reveals a mean annual vertical displacement of up to −5.3 cm/year, with limited horizontal motion during the study period from 2014 to 2019. 2.
Four large subsiding bowls are detected in four aquifer basins: Upper and Lower Kabul, Paghman, and Logar, with the former experiencing twice and three times as much subsidence as the three other sites. The maximum cumulative subsidence accounts for −23.8, −10.9, −12.9, and −8.6 cm, respectively.

3.
Analysis of groundwater level change illustrates its synchronized variations with subsidence: the larger the groundwater level drop, the greater the subsidence rates observed. Similarly, seasonal fluctuations in the groundwater level in Logar basin coincide with observed non-linear deformation patterns. Cross-wavelet transform analysis allows the estimation of the temporal lag of the system response (ground subsidence) to triggering factor (groundwater decrease) in Logar basin.

4.
A greater magnitude of subsidence is observed in areas with compressible clayey sediments in shallow geology.
The findings of the study suggest that ground deformation in the study area is largely anthropogenic rather than natural. In the context of climate change and increasing urban sprawl of Kabul, monitoring ground deformation using space-borne InSAR is an invaluable tool that provides new opportunities to inform future approaches for sustainable management of groundwater resources and improved planning of managed aquifer recharge operations. Within the framework of Terzaghi's consolidation theory, the synergistic use of InSAR-derived results with detailed hydrogeological and geotechnical data has proved itself a robust tool to provide a holistic picture of the mechanisms inducing the ground deformation and could potentially be applied over wide regions of interest. Wavelet transform analysis implemented with temporally-dense Sentinel-1 enabled a recurring clay shrink-swell hazard in the area from long-term subsidence phenomena to be distinguished.