Monitoring land subsidence induced by tectonic activity and groundwater extraction in the eastern Gediz River Basin (Türkiye) using Sentinel-1 observations

, and subsequent land subsidence. This study aims to evaluate the roles of tectonic activity and groundwater withdrawal in land-subsidence and investigate additional factors like faults and soft soil thickness. The P-SBAS algorithm was applied, using 98 and 123 Sentinel-1 SAR images in ascending and descending orbits, respectively, from 2016 to 2020. Independent Component Analysis separated long-term displacements from seasonal variations in the InSAR time series data. InSAR analysis showed displacement rates up to (cid:0) 6.40 cm/year. Results reveal a direct correlation between InSAR displacement and soft soil thickness, highlighting aquitard layer compaction due to groundwater withdrawal and piezometric head depletion as the primary causes of land subsidence. It was not possible to disentangle the multi-year displacements with different rates, i.e., to isolate the displacements caused by tectonic activity from displacements due to the soft soil compaction induced by groundwater withdrawal. Nonetheless, the analysis identified two spatiotemporal displacement trends: one linked to long-term, linear plastic compaction of the aquitard due to groundwater withdrawal and the other involving long-term displacements with seasonal rebounds caused by yearly water level fluctuations.


Introduction
The scarcity of groundwater reserves in some semi-arid regions, combined with economic and population growth have increased pressure on this valuable resource, leading to unsustainable exploitation rates in some cases.Consequently, modern societies have been forced to face numerous challenges related to ground deformation caused by the over-exploitation of groundwater (Galloway and Burbey, 2011).Land subsidence due to groundwater withdrawal is one of the most extensive types of geological phenomenon induced by human activity in the 20th century, affecting many major cities and regions in the world (Herrera-García et al., 2021).This widespread phenomenon is well-known and studied in China (Bian et al., 2014;Chen et al., 2021;Jiang et al., 2021;Shi et al., 2008;Zhou et al., 2017), central Mexico (Chaussard et al., 2014b;Cigna et al., 2021;Yalvac, 2020); Italy (Bonì et al., 2017;Poluzzi et al., 2019;Teatini et al., 2005;Tosi et al., 2009), and Mediterranean regions in Spain (Béjar-Pizarro et al., 2016;Bonì et al., 2015Bonì et al., , 2016;;Herrera et al., 2010;Tessitore et al., 2016).
Different techniques have been used conventionally for monitoring and measuring the subsidence, such as GNSS (benchmarks and continuous), levelling, extensometer, etc. (Chen et al., 2021).However, those methods present some limitations such as the low spatial density of measurement points which requires time consuming and expensive monitoring campaigns (Chen et al., 2021;Tomás et al., 2014).As an alternative, the differential interferometric synthetic aperture radar (DInSAR) technology has been developed and widely applied for earth observation with millimetre level precision and wide spatial coverage, producing results that are not affected by weather conditions (Sadeghi et al., 2021).This technology is used to obtain mean velocity maps and displacement time series.The Small BAseline Subset (SBAS) technique uses an interferogram network by linking multiple SAR images with small baselines (multi-master).Traditionally this technique considers distributed scatterers (DS) using a spatial average within a twodimension window, usually named multilook (ML), which allows for the reduction of noise at the expense of a resolution loss (Berardino et al., 2002;Minh et al., 2020;Sadeghi et al., 2021).
Türkiye is located on the Anatolian microplate, eastern Mediterranean, which is known as a very complex tectonic zone.The regional tectonic composition is an interaction between the Anatolian microplate and three other adjacent plates, namely the Eurasian, African and Arabian, which have led to numerous deformations and important tectonic activity, including earthquakes, particularly in the north Anatolian fault zone (Cavalié and Jónsson, 2013;Tatar et al., 2012).This tectonic activity has led to a horst and graben structure which dissects the Menderes massif into three units.These tectonic valleys, mostly located in the west of Türkiye and filled with detrital materials, are one of the regions most affected by land subsidence in the country (Dogan et al., 2022;Hacıoglu et al., 2021;Imamoglu et al., 2022).During the last two decades several studies have observed and monitored the faults which border the grabens, concluding that the normal faults induced structural damage via the differential subsidence along their traces (Imamoglu et al., 2022).However, other factors have been identified as potential cause of land subsidence, such as the overexploitation of the aquifers due to intensive agricultural and other anthropogenic activities (Orhan, 2021).The Alas ¸ehir-Sarıgöl sub-basin (ASSB) situated within The Gediz River Basin (GRB) is one of the abovementioned tectonic valleys; it is located in western Türkiye and developed from the regional extension horst-graben systems.Agriculture is the dominant economic activity in the area, accompanied by secondary economic activities such as animal husbandry, industrial, geothermal energy production, and mining; along with hosting densely populated municipalities (Duru et al., 2018).Since groundwater is the main source for agricultural and household water demand, the interaction between these water-intensive human activities and the natural features jeopardizes the quantity and quality of water resources, thereby rendering the Gediz basin as one of the most stressed basins in the country (Duru et al., 2018;Elçi et al., 2015).Regional tectonics also play an important role in the towns of the valley.The noteworthy reason is because some important towns are built over some active faults placed in the southern margin of the basin (e.g., Alas ¸ehir and Sarıgöl), where the graben presents more seismic activity and the highest earthquake record (Poyraz et al., 2019).The low-angle detachment fault limiting this margin has caused some structural damage in houses, streets and other infrastructures.Previous studies measured and monitored the subsidence rate by means of geodetic techniques such as InSAR and Global Navigation Satellite System (GNSS).Koca et al. (2011) calculated a cumulative vertical displacement in a Sarıgöl fault spot of − 60 to − 85 cm for the period 2000 to 2010.Poyraz et al. (2019) measured vertical displacements in different points within the valley using GNSS benchmarks during 2013 to 2015 and found a maximum displacement rate of − 90 mm/year along Sarigöl fault.Additionally, Poyraz and Hastaoglu (2020) monitored the valley by InSAR using TerraSAR-X SAR images acquired from May 2014 to November 2015, and observed a maximum displacement rate of − 50 mm/year in the Line of Sight (LOS) direction.The previous studies suggest that the source of the subsidence cannot be attributed to only the creep movement caused by tectonics, but that a component of the subsidence may be induced by groundwater extraction.
The aim of this study is to evaluate the role of tectonics and groundwater withdrawal in the Alas ¸ehir-Sarıgöl sub-basin on land subsidence and evolution of faults, as well as attempt to demonstrate that subsidence is linked to water extraction, and investigate other conditioning factors such as the soft soil thickness.To this aim, the parallel solution of the SBAS algorithm (P-SBAS), allocated in the Geohazard Exploitation Platform (GEP), was used to process 123 SAR images in descending orbit and 98 in ascending orbit acquired by Sentinel-1 between 2016 and 2021.Investigating land subsidence in the Alas ¸ehir-Sarıgöl sub-basin by combining P-SBAS and ICA provided important insights and understanding that can be extended and applied beyond this specific study area.The application of the ICA technique to InSAR data has rarely been applied in the past and the present study reveals the usefulness of this technique to identify subsidence trends and also to explain the main causes behind this phenomenon.Its implication involves worldwide significance and can contribute to the resolution of critical issues related to water resources, natural hazards and environmental sustainability, making it important for a global audience.This paper is organized as follows: Section 2 is dedicated to the geological and tectonic settings of the study area: Section 3 describes SAR images processing and the Independent Components Analysis (ICA) applied to the time-series; Section 4 shows the P-SBAS results and validation and the ICA results; Section 5 discusses the effect of different factors, such as fault locations, soft soil thickness, seasonality and changes in the piezometric level on land subsidence.Finally, the main conclusions are presented in Section 6.

Study area
The study area is located in the ESE-WNW trending Gediz graben, also known as the Alas ¸ehir Graben, and encompasses sections of the Gediz River and the Alas ¸ehir Creek, a tributary of the Gediz River (Hacıoglu et al., 2021;Fatih Poyraz and Hastaoglu, 2020;Üner and Dogan, 2021) (Fig. 1a).The ASSB exhibits a typical Mediterranean climate, with hot and dry summers and cold and rainy winters; the average temperature is about 15.2 • C with an annual precipitation of approximately 617 mm (Elçi et al., 2015).Agriculture dominates the land use in this part of the basin covering more than 80% of the area, with vineyards representing the maincrop type occupying approximately 70% of the cultivated surface.The main residential areas are concentrated in Salihli, Alas ¸ehir, Yes ¸ilyurt and Sarıgöl, with an approximate population of more than 250,000 inhabitants.The continuous urban fabric corresponds to approximately 1.8 km 2 and discontinuous urban fabric to 1542 km 2 (Copernicus, 2018).

Geological setting
The Alas ¸ehir-Sarıgöl sub-basin placed in the Gediz River Basin alluvial plane develops on a horst-graben system, the graben is filled with roughly 500 m of Pliocene and Quaternary alluvial material underlain by more than 2 km thick Miocene alluvial and lacustrine sedimentary rocks (Sözbilir, 2002;Üner and Dogan, 2021).This is one of the youngest and the most developed continental extensional basins in western Anatolia (Poyraz et al., 2019).The rock units in this graben are divided into two main groups: (a) metamorphic rocks that correspond to the basement and outcrop on the borders of the valley.They consist of Precambric-Palaeozoic Menderes Massif, composed by mica schist, phyllite, meta-sandstones and crystallizes limestone (Haluk Selim and Yanik, 2009); (b) sedimentary cover, the contact between the basement rocks and the sedimentary fill is defined by normal faults in the margin and unconformably underneath the graben (Hacıoglu et al., 2021;Üner and Dogan, 2021).Fig. 1b shows in a light-yellow colour the Quaternary material filling the basin (uppermost unit).This material is mainly composed of fluvial sediments of the modern Gediz River.Meanwhile, Neogene basin fill (Miocene and Pliocene) corresponds to sedimentary sequences, which are composed by fluvial-alluvial, alluvial-fan, debris and clay rich interval layers with very low permeability.The Alas ¸ehir formation corresponds to the lowermost and oldest layer of the graben fill and it is composed by conglomerate, sandstones, siltstone and shales with high organic material content (Hacıoglu et al., 2021;Üner and Dogan, 2021).

Tectonic and seismic activity in the Gediz Basin
The Anatolian plate is one of the most seismically active regions in the world, due to the interaction with the Eurasian, African and Arabian plates (Poyraz et al., 2019).Western Anatolia is located in the rapid continental crustal Aegean extension region (Haluk Selim and Yanik, 2009) developed in a back-arc tectonic environment, which is undergoing a roughly north-south extensional deformation, with velocity values between 1.5 and 3.0 cm/year (Hacıoglu et al., 2021;Poyraz and Hastaoglu, 2020) as a result of north-dipping Hellenic subduction zone in Africa-Eurasia convergent boundary (Ozdemir et al., 2021).This region has been affected by several destructive earthquakes, at least 13 earthquakes occurred in the past, with most of them originating in the Gediz graben causing important damage in residential areas (Poyraz andHastaoglu, 2020) (e.g., Salihli earthquake in 1965, Mw 5.8;Alas ¸ehir earthquake in 1969, Mw 6.5, and Gediz earthquake in 1970, Mw 7.2 (Haluk Selim and Yanik, 2009).The Gediz graben is part of the Aegean horst-graben system (Haluk Selim and Yanik, 2009), limited by the currently inactive low-angle (~10 • ) Gediz detachment fault to the south, which is the most important structural element of the system, separating the Neogene sediments from the metamorphic basement corresponding to Menderes Massif (Hacıoglu et al., 2021;Haluk Selim and Yanik, 2009;Poyraz et al., 2019).A second structural and important element in the study area are the E-W oriented high angle oblique and active faults that cross the detachment faults (Sözbilir, 2002;Üner and Dogan, 2021) as a result of horst-graben system development.Currently, most seismic activity occurs along these faults (Poyraz and Hastaoglu, 2020).

Methods and data processing
The methodology applied in this work is summarized in Fig. 2 and can be divided into two main sections.Section 1 is related to the InSAR processing, calculation of 2D components and validation of the accuracy of the displacement measured by SAR image acquisitions using a Continuous GNSS (CGNSS).The velocity map obtained by the P-SBAS algorithm is used for the geological interpretation and the analysis of the relationship of the subsidence with the soft soil thickness, and the active faults within the basin.Section 2 corresponds to the analysis of the time series derived from the InSAR processing by the application of Independent Component Analysis (ICA).The objective was to identify the different deformation trends such as seasonal behaviour, and their potential relation with factors such as precipitation and/or the aquifer recharge and discharge, as well as to determine linear deformations that could be related to the permanent deformations of the aquifer associated with lowering of the piezometric level.

P-SBAS
Differential InSAR (DInSAR) is a technique that measures the ground deformation generating interferograms using two complex SAR images (Ferretti et al., 2007), but incorporating an external Digital Elevation Model (DEM) in order to subtract the topographic contribution from the interferometric phase, described in Eq. ( 1), where Δφ is the total interferometric phase, Δφ flat is the flat-earth component due to range distance difference between pixels, Δφ def is the relative displacement,Δφ topo is the topographic component,Δφ atm is a phase component related to the atmospheric phase screen and Δφ noise comprises the errors due to different decorrelation sources (Declercq et al., 2017;Tomás et al., 2010;Zebker and Villasenor, 1992).
Nevertheless, DInSAR is a method which does not resolve other limitations related to the interferogram phase, such as temporal decorrelation, spatial decorrelation and atmospheric delay (Gao et al., 2016;Hooper, 2008).Thus, to address all this issues more advanced algorithms have been developed during the last two decades, such as, the Permanent Scatters (PS) and Small BAseline Subset (SBAS) techniques (Casu et al., 2006;Hooper, 2008;Minh et al., 2020;Yalvac, 2020).
SBAS is an Advanced DInSAR technique (A-DInSAR), which has been developed to minimize the temporal and geometric decorrelation and noise effects.It enables the production of deformation velocity maps and displacement time series with centimetre to millimetre accuracy (De Luca et al., 2018;Yalvac, 2020).In order to reduce the noise, this technique uses a multi-reference interferogram network characterized by short temporal and perpendicular baselines (Berardino et al., 2002;Manunta et al., 2019).The P-SBAS phase unwrapping procedure guarantees fully consistent solutions when dealing with multilook interferograms sequences characterized by short-time temporal baseline (De Luca et al., 2021) The use of small (temporal) baseline interferograms ensures high coherence between the pairs, since extensive changes on the surface due to land cover or human activities (e.g., crop growth) could cause temporal decorrelation (Nefeslioglu et al., 2021).
An improved version of the SBAS algorithm, called Parallel SBAS (P-SBAS), was developed to better exploit the distributed high performance computing infrastructure in an efficient way.P-SBAS allows the use of multi-core and multi-node programming techniques, guaranteeing supportable scalable performance in order to process massive amounts of data (De Luca et al., 2018;Zinno et al., 2017).This approach was run on the Geohazard Exploitation Platform (GEP), which is an on-demand web tool launched in 2015 by the European Space Agency (ESA) (Bru et al., 2021;Cigna et al., 2021;Ruiz-Armenteros et al., 2018).This platform can be used to perform DInSAR analyses with a friendly interface, enabling the users to run P-SBAS algorithm with the aim of monitoring different phenomena, such as earthquakes, landslides, land subsidence, and volcanic activity (Ikuemonisan et al., 2020).
In the ASSB, the P-SBAS algorithm was performed on 123 Sentinel-1 SAR images in descending orbit (from June 2016 to July 2020) and 98 in ascending orbit (from December 2017 to February 2021).For the P-SBAS processing the coherence threshold was set at 0.80 for descending and ascending modes.A threshold of stability of ±0.7 cm/year was established for the P-SBAS processing.Land cover in the study area mainly includes: urban fabric in cities and towns that exhibit high coherence among the SAR scenes; agricultural crops with low coherence; and mountain areas with scarce vegetation and high coherence.
The reference point is located in the deformation far field, more precisely within the mountainous region characterized by gneiss and clastic sedimentary rocks, which surrounds the alluvial valley and is in close to the stable GNSS point "KZLC."This station serves as a benchmark GNSS with only three available measurements from 1995 to 2016, during which the cumulative deformation reached 14 mm.This justifies the choice to consider this location as the reference point.It is crucial to emphasize that the correct selection of the reference point is paramount for the analysis.This point functions as a fixed anchor against which we can accurately measure the relative motion of other points in the study area.This stable reference point establishes a dependable baseline for interferometric measurements, aiding in the identification of deformation patterns, such as land subsidence or uplift.The absence of a stable reference point can make it challenging to differentiate between actual ground motion and potential data errors or noise.Consequently, choosing an unstable reference point can result in inaccurate interpretations of subsidence or deformation rates.Moreover, as per the P-SBAS manual, the algorithm automatically refines the reference point and selects the one with the best coherence conditions near the userselected point.Additionally, at the end of the processing, the algorithm implements an average reference across the entire scene to mitigate dependence on a single point.This approach helps alleviate the impact of atmospheric noise on the reference point, effectively creating a form of 'absolute' reference for all zero-mean points.

2D component calculation
The combination of ascending and descending satellite datasets of ASSB enables calculation of the 2D component assuming that SAR is not sensitive to north-south horizontal motion (i.e., V n = 0).Hence, combining the ascending and descending LOS vectors is possible to carry out the geometry decomposition using Eqs.( 2) and (3) (Béjar-Pizarro et al., 2017;Cigna et al., 2019).Where: H d , H a are the vertical directional cosine of descending and ascending LOS unit vector, respectively.E d , E a are the east-west directional cosine of descending and ascending LOS unit vector, respectively, and vLOS d , vLOS a are the mean velocity in the Line of Sight (LOS) of the descending and ascending satellite orbit passes, respectively.
For this purpose, firstly, we calculated the common velocity period for both datasets (ascending and descending) between December 2017 to July 2020.Secondly, a 100 m buffer around each measuring point (MP) of ascending orbit was built and used for the selection of the common measurement points with the descending orbit.A new dataset was produced containing only the common MP for both orbits, and used for the construction of six different raster maps, vLOS a ; vLOS d ; H a ; H d ; E a ; E d using the inverse distance weighted (IDW) method with a 100 m searching radius and minimum of two data points, in order to calculate the eastward and vertical component using raster calculation GIS tool and Eqs. ( 2) and (3).
Since this methodology is suitable only for the common MP in both orbits the density of MPs decreases considerably, and an important amount of information is lost.This has presented a challenge in this study, considering the fact that the ASSB is predominately an agricultural area with a low number of persistent scatterers when processed with a high coherence threshold.Consequently, the number of MPs was increased within the valley by including additional MPs from one of the available datasets originally rejected during the calculation of the 2D displacement components.The displacements of these MPs were then projected to the vertical direction neglecting the horizontal movement.It should be noted that this hypothesis is valid for most of the flat subsiding areas in which horizontal strains associated to land subsidence are not significant (rarely as large as 2 ppm) (Galloway and Burbey, 2011).

Principal and independent components analysis
Interferograms can be considered as a mixture of signals (Eq. ( 1)).In this study Principal Component Analysis (PCA) and Independent Component Analysis (ICA) have been applied with the purpose of separating and identifying spatiotemporal patterns of long-term deformation and seasonal variations.PCA is a well-known statistical method used to analyse time series data aimed at finding the linear transformation (principal component, PC) which best fits the maximum variance from observations (M) and variables (N) (Chaussard et al., 2014a;Cohen-Waeber et al., 2018;Karimzadeh and Matsuoka, 2020), to achieve a reduction in the dimensionality of a dataset, while at the same time preserving as much variability as possible (Jollife and Cadima, 2016).The displacements measured by P-SBAS technique were used as input for the PCA analysis, considering only the time-series of the MPs located within the sub-basin for the descending orbit (~14,000 MP).The time series of the selected MPs (observations) were used to build a matrix, which was used as input in the general equation of PCA (Jollife and Cadima, 2016;Karimzadeh and Matsuoka, 2020): where, X is the data matrix including the time series, T is the matrix of PC scores, P′ is the transposed matrix with the loading factors, and r is the residual matrix (Jollife and Cadima, 2016;Karimzadeh and Matsuoka, 2020).
The major advantage of PCA is to order the uncorrelated variables, thereby capturing the variance in a more efficient way.However, this method may miss key information by the mixture of some part of the data trend (Cohen-Waeber et al., 2018).In addition, PCA may fail for large datasets, forcing the users to reduce the number of pixels analysed (Gaddes et al., 2018;Karimzadeh and Matsuoka, 2020).Meanwhile, ICA is a statistical computational signal technique which assists with separating linearly mixed independent sources, rather than extracting components based on signal variance (such as done with PCA), assuming that the components do not follow Gaussian distributions (Chaussard et al., 2017;Ebmeier, 2016;Wang et al., 2021).Similar to PCA, ICA tries to find the independent sources, but from independence point of view, ICA is a more robust processing than the uncorrelatedness searched by PCA (Gaddes et al., 2018).Since the amount and the order of ICs to maintain is given by the interpretation of the explained variance obtained from the PCA, ICA can be considered as an extension of PCA (Estelle Chaussard et al., 2017;Estelle Chaussard and Farr, 2019).Furthermore, several ICA algorithms including FastICA (Hyvärinen, 1999) which was used in this study, combine both approaches.FastICA applies PCA in the first step in order to isolate the ICs and also for data whitening, which is usually required and incorporated into the dimension reduction step.Data whitening involves the transformation of the original data to enhance its suitability for Independent Component Analysis (ICA).This process comprises three key steps: centring, decorrelation, and scaling.The primary objective of data whitening is to establish uncorrelated data with unit variance.When applied in the FastICA algorithm, this pre-processing step simplifies the search for statistically independent components, thereby enhancing the effectiveness of the algorithm in extracting meaningful information from the original data (Ebmeier, 2016;Gaddes et al., 2018).The main results of the FastICA analysis are the eigenvalues, which correspond to the diagonal values of the correlation matrix; the eigenvectors, which are the magnitude of the contribution of each original variable to each PC or IC; and the IC score maps, which illustrate the spatial distribution of each component and correspond to the contribution of the original matrix at every point (Chaussard et al., 2014a;Cohen-Waeber et al., 2018).

P-SBAS results
The ground velocity deformation maps derived from P-SBAS processing in the study area are displayed in Fig. 3. Descending and ascending orbit results are depicted in Fig. 3a and b, where the red points (i.e., the negative values) are the areas moving away from the sensor along the LOS and the blue points (i.e., the positive values) represent the areas moving toward the sensor.The P-SBAS results cover an area of around 13,000 km 2 for the processed descending images and 24,000 km 2 for the ascending.A total of 1,047,952 and 1,689,292 MPs were selected for descending and ascending orbit, respectively, resulting in a pixel density of 81 and 70 pixel/km 2 , respectively.Peak LOS rates within the study area have a similar range, varying from − 6.37 to 1.84 cm/year for the descending orbit and from − 6.40 to 2.87 cm/year for the ascending.From Fig. 3 the low MP density within the sub-basin can be observed.This is due to the type of land use; this area is strongly affected by quick changes caused by agricultural activity and crop phenology over time, reducing the coherence and hence the spatial density of points.The results of both datasets identify the same areas of maximum subsidence, which are concentrated in the centre of the basin and on the southeast edge of the tectonically active graben (Sarıgöl town).Land movements toward the satellite are concentrated on the northwest edge of the basin (black square).The comparison of the uplifting spot with the geological information suggests that this positive deformation might be interpreted as a local extension dynamic (newest horst), since this small area is composed of Miocene sedimentary rock surrounded by younger Pliocene rocks.Fig. 3c shows the vertical displacement calculated following the methodology described in Section 3.1.2,the results exhibit a subsidence distribution similar to LOS rate maps, where the main vertical displacement are located within the basin and concentrated in the southeast border.Meanwhile, the horizontal displacements shown in Fig. 3d capture the regional tectonic trend movements toward west direction (points in red and orange) related to the continental crustal Aegean extension.These MPs are identified not only in the sub-basin, but also in the mountains.Points in blue are representative of east horizontal movement, with all of these MPs located within the graben.The main of these displacements occurs along the dip angle faults, for that reason they might be related more to a local tectonic linked with the Gediz graben evolution.

P-SBAS validation
To validate the results from satellite data, the time series of an available in-situ CGNSS station (SALH) and nearby InSAR MPs were compared (Fig. 3a and b).This station is located in the city of Salihli, covering the time period between November 2008 to October 2020.
InSAR dataset validation was achieved using the ValInSAR code (Navarro-Hernández et al., 2022) comparing the available data from the Continuous GNSS station with the nearest MP.It should be noted that only the descending dataset was compared with the GNSS time series, since the ascending dataset did not include any MPs in the area close to SALH station (Fig. 3b).In order to compare the InSAR dataset, a 150 m diameter buffer was delineated around the GNSS location that defines a zone for averaging the InSAR data pixels.Then, the GNSS and InSAR time series were compared for the common period as shown in Fig. 4a.Fig. 4b displays the correlation plot between Sentinel-1 and CGNSS acquisitions.The statistical results of the comparison between CGNSS and InSAR (P-SBAS descending orbit) time-series are also included in Fig. 4. For the comparison between Sentinel-1 acquisitions and SALH station the RMSE is 0.91 cm, the mean error is 0.76 cm, the standard deviation of the error is 0.50 cm, the maximum discrepancy is 2.33 cm, and the coefficient R 2 value is 0.91 (signifying good forecasting).These results demonstrate a good agreement between both sources of data for SALH station since both techniques show a similar subsidence trend over time.
Recognizing the need for a minimum of two ground control points (GCP) to improve validation accuracy, and due to the inability to compare InSAR results with in-situ data in our study area, we adopted the methodology proposed by Crosetto et al. (2021) for the validation of the Ground Motion Service (EGMS) datasets.This approach involved specific criteria:  Crosetto et al. (2021) suggest that for land subsidence three kind of validation activity (VA) can be carry out: VA3 "Comparison with a ground motion services"; VA4 "Comparison with inventories of phenomena, event, damage or anthropic activities" and VA7 "Comparison with in-situ monitoring data, including ground or single building motion monitoring or groundwater level measurement".Since no ground motion service exists within the ASSB, only two validation activities were applied for this pilot site.
To further validate the findings, the authors conducted a comparison of InSAR results with GNSS benchmarks (specifically, BGCL and AZKL points) from a prior study (Poyraz et al., 2019) that are located within the sub-basin but monitored during a different time period (2013)(2014)(2015).Despite not being measured during the same time period, the benchmarks exhibited only minor differences in vertical velocities, with BGCL at − 0.67 cm/year and AKZL at − 0.57 cm/year.
Lack of in-situ measurements, only additional information of the basin can be compared with the spatial distribution and temporal evolution of the displacements to estimate their consistency.In the Alas ¸ehir-Sarıgöl sub-basin the available information is related to the mapped active faults, the distribution of the soft soils thickness, and the detection of damages in the urban areas.The results of this analysis will be detailed in Section 5.

PCA/ICA results
Applying FastICA algorithm to the P-SBAS results made it possible to determine the first order trends obtained by three uncorrelated independent components (ICs), which accounted for more than 99% of the eigenvalues (Fig. 5a).The ICs correspond to a transformation of the time-dependent information into a dataset of representative spatial pattern of deformation, which is displayed as score maps; these score maps emphasize a signal coherence in space, but not over time (Fig. 5c-5d).The distribution of IC scores provides an understanding of the correlation between the trend of each IC and the measuring point (MP).IC1 explains nearly 96% of the variance (Fig. 5a) and shows a spatial signal distribution (Fig. 5c) similar to the mean velocity map (Fig. 3a), which suggests that this component corresponds to the long-term deformation (MPs in red colour on Fig. 5c).IC2 explains 2.5% of variance, showing a spatial signal in the southeast and centre of the basin (Fig. 5d).This component recovers the seasonal deformation behaviour shown in Fig. 5b.As it can be seen, the peaks are in agreement with the seasonal cumulative rainfall.However, despite the seasonality related to the precipitation, there is a deformation component which continually decreases over time.This component might be correlated to the piezometric head changes in the deeper aquifer, and will be discussed in Section 5.2.IC3 explains approximately 1% of the variance and shows a spatial extensive deformation signal for all of the basin, recovering a seasonal deformation pattern (Fig. 5b) with trough and peaks in agreement with the changes in seasonal cumulative rainfall values.This trend probably encompasses the changes produced in the soil moisture of the surface, which is mostly composed of granular material.This kind of soil material typically responds faster to the loading and unloading cycles.This is supported by Fig. 5e, where it can be seen that the majority of the valley surface shows a IC3 positive score and is in agreement with the location of the Quaternary alluvial material displayed in the geological map.

Discussion
The investigation of land subsidence in the ASSB using advanced techniques like P-SBAS and ICA with the purpose of identifying the roles of tectonic activity and groundwater withdrawal can be of global significance, as these issues are not isolated but are actually a common problem in many regions worldwide which are facing similar challenges (e.g., cities in central Mexico with horst and graben valleys dealing with groundwater exploitation problems (Castellazzi et al., 2016;Figueroamiranda et al., 2018)).Thus, understanding the interaction of both factors and the impact in the ASSB can provide valuable insights applicable to similar regions globally.Furthermore, the global trend of rapid urbanization often relies on groundwater.Therefore, understanding the link between groundwater extraction and subsidence is vital for urban planning and infrastructure development worldwide.Land subsidence also impacts ecosystems, agriculture, and water quality, making it an international concern for environmental conservation (Bonì et al., 2022).
Policymakers and decision-makers worldwide require reliable data to formulate effective policies, making the ASSB research crucial in informing decisions related to groundwater management and territorial planning (Bonì et al., 2022).Additionally, given the substantial economic impact of subsidence on infrastructure, agriculture, and resource management, understanding its drivers in the ASSB can guide economic decisions in regions facing similar challenges.
The following subsections discuss the key findings regarding the relationships between land subsidence and the various factors analysed in this study.

Relationship with soft soil thickness
Displacements (δ) are generally evaluated in geotechnical engineering by considering three different parameters (Tomás et al., 2010), first, the thickness of the potentially deformable soil (D); second, the variation of the effective stress state, which in this case can be expressed as an effective piezometric level drop (Δh); and a soil modulus (S), which depends on the stress state, and represents the deformability of the soil.This relationship can be expressed as follows: According to this expression, a direct relationship exists between subsidence and soft soil thickness, such that the thicker the accumulation of soft soil (D), the larger the displacement.The bedrock structure can also have a strong influence on the soft soil thickness (D).Buried faults associated with grabens, horsts and other geological structures can cause sudden spatial changes in the soft soil deposits.
Fig. 6a displays the soft soil thickness map for the ASSB.This map was developed using the lithological descriptions from 118 boreholes mostly located in the alluvial zone.As the aquifers with high groundwater yield occur in these unconsolidated units, the mechanical response to the increase in the effectives stress is one of the components which controls the evolution and the magnitude of the displacements.The distribution of the soft soil thickness is shown in Fig. 6a.There is evidence of some agreements between high soft soil thickness and the highest displacement areas detected by InSAR located to the east and in the centre of the valley (red circle in Fig. 6a and b).To analyse the relationship between soft soil thickness and subsidence, a spatial correlation was investigated by comparing the geological information of the boreholes with the vertical displacement rate information obtained from P-SBAS.To this aim, all the MPs within a 150 m buffer around the boreholes were considered for this comparison to obtain the correlation (Fig. 6c).The correlation plot exhibits a good agreement between both variables with a R 2 of 0.853 despite of the small number of samples due to the reduced number of MPs located inside the sub-basin.This correlation shows that the highest subsidence rates (displacement rate of up to − 4 cm/year) develop in the areas of maximum soft soil thickness (from 50 to 100 m) and vice versa.These results demonstrate that the thickness of the graben fill sediments play a significant role in the occurrence of land subsidence and can be considered as a potential predisposing factor in the ASSB.
Fig. 7a displays the IC1 score map (Section 4.3), IC1 explains the 96% of the eigenvectors and its signal has a similar spatial distribution as the LOS deformation rate map.The positive IC1 score values can be correlated to the increase of clay material thickness, mainly accumulated in the centre of the basin.Similarly, the spatial distribution of PC1 can also be correlated with the soft soil thickness.Additionally, the comparison between P-SBAS time series and eigenvector time series within zones 1 and 2 (Fig. 7b and c) reveal that the higher positive IC1 score values are matching with the higher subsidence areas.This suggests that this PC corresponds to the main trend of deformation caused by long-term land subsidence.Furthermore, this component is strongly related to the compaction of the soft soil material due to the increasing effective stress.The over-exploitation of groundwater and resulting continuous decrease of piezometric levels would explain the increase in the effective stress.

Relationship with seasonality and groundwater extraction
Regarding the use of groundwater, approximately 86% of the total groundwater withdrawals are used for irrigation of the agricultural areas, while the remaining 14% is allocated to urban water demand and husbandry activities.The fraction of groundwater withdrawal for industrial purposes is relatively small (GDWM, 2018).Eq. ( 5) suggests that the effective piezometric level drop (Δh), which causes soil effective stress increase, has a direct relationship with subsidence.Therefore, the higher the piezometric level drop, the higher the ground surface deformation.However, we should consider that Eq. ( 5) is valid for a complete dissipation of excess pore pressure, and the dissipation of excess pressure will extend in time depending on the soil hydraulic conductivity.Since piezometric level drop (i.e., Δh) is one of the responsible factors for land subsidence, a map displaying the change in piezometric level was obtained using available observation data from 11 piezometers (Fig. 8c), by considering that the distribution of pumping wells strongly conditions the groundwater drawdown.This map suggests that the main piezometric level drop occurs in the centre of the plain due to groundwater withdrawal for irrigation and also along the southern border where the largest towns in the area are located (i.e.,    Salihli, Alas ¸ehir and Sarıgöl).Comparing this information with the IC2 score map in Fig. 8d it can be seen that they are in agreement and that the distribution of IC2 overlaps in some areas with the IC1 extent, suggesting that aquitard layer compaction can be related to the groundwater extraction.Additionally, in order to analyse the correlation with groundwater withdrawal and IC2, a comparison between the IC2 eigenvector and the water evolution was carried out as can be seen in Figs.8a-8b.This comparison shows a good agreement between the IC2 eigenvector time series and hydraulic head changes for the wells 55,129 and 22,043 from 2016 to 2018.It should be noted that both time series reveal no lags (or very short delays, but it is difficult to confirm due to the sporadic groundwater level data).Moreover, the IC2 time series displays a long-term subsidence with significant seasonal rebound, showing uplift and subsidence period also correlated with the cumulative seasonal rainfall peaks in Fig. 5b but with monthly delays between 2016 and 2018.

Relationship with active faults
Unfortunately, it was not possible in this study to separate the multiyear deformation with different rates by the application of ICA, i.e., to isolate the deformation caused by tectonic activity from deformation due to the soft soil compaction induced by groundwater withdrawal.Differential ground subsidence has been defined before in previous studies as a gradual settlement triggered by groundwater withdrawal, which leads to discontinuities and ground collapses according to the direction of a controlling tectonic structure (Figueroa-miranda et al., 2018).These aligned structures have been detected in the study area by the Sentinel-1 satellite and they agree with the mapped active faults, as can be seen in Sarıgöl fault (Fig. 9a).In Fig. 3a it can be noticed that the graben is more active at the southern edge.Previous researchers detected significant seismic activity along the southern margin, where the dip-slip normal faults in E-W direction are serving as a contact between the Miocene, Pliocene and Quaternary rocks and the metamorphic basement rocks (Poyraz et al., 2019).
Usually, graben structures control the thickness of the sediments accumulated in the basin, and therefore a greater number of unconsolidated materials has accumulated in the hanging wall of the fault compared to the filled material accumulated on the footwall.This effect is observed in Fig. 9b which plots different cross-sections (A-A' to D-D′) across the fault, showing a sharp change in the velocity of the deformation on both sides of the Sarıgöl fault.In this area the thickness of the clay and silty material changes drastically on either side of the fault due to the topography of the bedrock causing a differential compaction and, consequently resulted in major structural damage in the city (Figs.9c-d).
It should be also noted that the bedrock faults can also abruptly disrupt the lateral continuity of the aquifer layers, influencing the differential compaction at both sides of the fault.

Conclusion
In the Alas ¸ehir-Sarıgöl sub-basin of the GRB (Türkiye), the P-SBAS processing applied to Sentinel-1 SAR images between 2016 and 2020 has revealed that the maximum subsidence rates are mainly concentrated in agricultural and urban areas.InSAR results also exhibit a direct relationship with soft soil thickness: the thicker the soft soil thickness, the higher the subsidence rate.Additionally, there is a relationship between the piezometric decline and the subsidence in those areas in which high groundwater depletion coincides with areas of high soft soil thickness.ICA results reveal two types of spatiotemporal deformation trends captured by the two first components: a) IC1, which corresponds to long term and quasi-linear deformation due to the compaction of the aquitard, and b) IC2 which represents the long-term deformations with seasonal rebounds produced by the seasonal loading and unloading cycles due to water level fluctuations.In a general sense, the results found in this study indicate a human-intensified compaction of the basin filling due to groundwater withdrawal and piezometric head depletion strongly controlled by the graben structure.In fact, important differential compaction rates have been measured on both sides of the faults placed in the south edge of the basin inducing earth fissures (e.g., in Sarıgöl town) and causing important structural damage to the urban areas.
The P-SBAS processing applied to Sentinel-1 SAR images between 2016 and 2020 has revealed an important subsidence case in Alas ¸ehir-Sarıgöl sub-basin (Türkiye), with maximum rates of more than 5 cm/ year, affecting both urban and agricultural land uses.From our study, we can conclude a direct relationship between the InSAR results and soft soil thickness.Additionally, a clear relationship was also identified with the piezometric pressure drop in those areas, where high groundwater pumping rates coincide with areas of thick soft soil.
The novel application of the ICA to InSAR data let us to identify two types of spatiotemporal deformation trends captured by the two first components: a) IC1, which corresponds to long term and quasi-linear deformation due to the compaction of the aquitard, and b) IC2 which represents the long-term deformations with seasonal rebounds produced by the seasonal loading and unloading cycles due to water level fluctuations.These two main components were able to explain almost 98% of the eigenvalues.
The results indicate a human-intensified compaction of the basin filling due to groundwater withdrawal and piezometric head depletion strongly controlled by the graben structure.In fact, important differential compaction rates have been measured on both sides of the faults placed in the south edge of the basin inducing earth fissures (e.g., in Sarıgöl city) and causing important structural damage to the urban areas.
Leaders and decision-makers around the world depend on credible data to formulate effective policies.This research carries significant importance for making decisions related to groundwater management and land-use planning.Additionally, given the impact of land subsidence on infrastructure, agriculture, and resource management, understanding its underlying causes in the Alas ¸ehir-Sarıgöl sub-basin can provide valuable insights for decision-makers in regions facing similar challenges.In the future, the fundings derived from this research, can be integrated to subsidence modeling and it will enable authorities to predict land subsidence under different scenarios and supporting to formulate corresponding mitigation strategies.

Fig. 1 .
Fig. 1.(a) Location of the ASSB (red polygon).(b) Geological map of the study area with delineation of regional faults.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Flowchart of the methodology applied in this work.

Fig. 3 .
Fig. 3. P-SBAS displacement rate and SALH CGNSS location.(a) LOS displacement rates for descending orbit, (b) LOS displacement rates for ascending orbit, (c) vertical displacement rates, V U; and (d) east-west horizontal displacement rates, V E .
a. Completeness and consistency of the InSAR products: This criterion involved qualitatively evaluating the consistency of the InSAR results with the environmental, urban, and geological context (e.g., landforms and geomorphology) of the validation site.b.Accuracy of the InSAR products: This criterion assessed the ability of the processed InSAR products to closely approximate the true measurements, represented by higher-quality data obtained from in situ instruments or other remote sensing sources.It indicated that the processed InSAR products could accurately detect active motions and track their temporal evolution.

Fig. 4 .
Fig. 4. Validation of the descending orbit P-SBAS dataset: (a) comparison between SALH GNSS and P-SBAS time series.(b) Correlation plot between SALH CGNSS and P-SBAS time series data.Continued lines represent the 1:1 line, and dotted lines parallel to the 1:1 line represent the ±1 cm/error.

Fig. 5 .
Fig. 5. FastICA results: (a) Variance explained by each component in PCA, (b) ICs eigenvector time series compared with seasonal cumulative rainfall (blue), (c) IC1 score map, (d) IC2 score map, (e) IC3 score map.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.Relationship between soft soil thickness and land subsidence: (a) Soft soil thickness map derived from boreholes.(b) Vertical component of LOS displacement rate.(c) Correlation plot between soft soil thickness and vertical displacement rate.

Fig. 7 .
Fig. 7. Comparison between IC1 and LOS displacement at two locations with high soft soil thickness and land displacements: (a) time series comparison at zone 1.(b) Time series comparison at zone 2.

Fig. 8 .
Fig. 8. Relationship between seasonality and piezometric head evolution: (a) IC2 and piezometric head time series at the well 55,129.(b) IC2 and piezometric head time series at the well 22,043, (c) spatial distribution of groundwater level evolution from 2016 to 2018 and (d) IC2 map score.

Fig. 9 .
Fig. 9. Relationship between active faults and land subsidence: (a) Rate displacement map covering Sarıgöl fault influence area.(b) Analysis of vertical displacement sections A-A', B-B′, C-C′ and D-D′ along the Sarıgöl fault trace, (c), (d) and (e) Structural damages caused by the displacements in the town of Sarıgöl.