Use of InSAR data for measuring land subsidence induced by groundwater withdrawal and climate change in Ardabil Plain, Iran

The Ardabil plain, with an approximate area of 1097.2 km2 in northwestern Iran, has experienced land subsidence due to intensive groundwater withdrawal and long seasons of drought in recent years. Different techniques have been used to investigate and evaluate subsidence in this region including: Global Positioning Systems (GPS), Levelling, and Geotechnical methods. These methods are typically expensive, time-consuming, and identify only a small fraction of the areas prone to subsidence. This study employs an Interferometric Synthetic Aperture Radar (InSAR) technique to measure the long-term subsidence of the plain. An open-source SAR interferometry time series analysis package, LiCSBAS, that integrates with the automated Sentinel-1 InSAR processor (COMET-LiCSAR) is used to analyze Sentinel-1 satellite images from October 2014 to January 2021. Processing of Sentinel-1 images shows that the Ardabil plain has been facing rapid subsidence due to groundwater pumping and reduced rainfall, especially between May 2018 to January 2019. The maximum subsidence rate was 45 mm/yr, measured at the southeastern part of the plain. While providing significant advantages (less processing time and disk space) over other InSAR processing packages, implementation of the LiCSBAS processing package and its accuracy for land subsidence measurements at different scales needs further evaluation. This study provides a procedure for evaluating its efficiency and accuracy for land subsidence measurements by comparing its measurements with the results of the GMTSAR and geotechnical numerical modeling. The results of geotechnical numerical modeling showed land subsidence with an average annual rate of 38 mm between 2006 and 2020, which was close to measurements using the InSAR technique. Comparison of the subsidence measurements of the Ardabil plain using the LiCSBAS package with results obtained from other techniques shows that LiCSBAS is able to accurately detect land deformation at large scales (~ km). However, they may not be optimized for more local deformations such as infrastructure monitoring.


Description of study area and research methodology
Description of study area. The Ardabil plain, with an approximate area of 1097.23 km 2 , is an inter-mountain plain located in northwestern Iran with geographical coordinates of 38°5' to 38°27' N and 48°9' to 48°37' E ( Fig. 1). The Ardabil plain has been formed out of Quaternary alluvial deposits due to weathering and erosion of surrounding mountains 48 . Conglomerate, volcanic ashes, and lahars are the predominant lithological outcrops in the west of the plain. There are abundant springs in these units which infiltrate into the groundwater and affect aquifer recharge. In the past, most springs, qanats, and rivers were used for water supply, but now water supply is mostly based on the exploitation of wells. Based on the results of geophysical studies, pumping test data, and drilling logs, Ardebil aquifer has a maximum thickness of 220 m and is mainly composed of gravel, sand, and a small amount of clay. The transmissivity of the aquifer varies between 50 and 2200 m 2 /day, and the specific yield ranges from 0.021 to 0.14. Storativity along with transmissivity and specific yield are key properties of the groundwater system. They are typically obtained in situ from pumping tests and reflect the response of aquifers and aquitards to groundwater head changes. Knowledge of their values is necessary to quantify available groundwater resources 49,50 .
The general groundwater flows from other directions into the north-west of the plain 48 . More than 85% of exploited groundwater is used for agricultural activities, 14% is used for drinking purposes and 1% for industrial use. The yearly maximum and minimum groundwater levels are usually recorded in May and September, respectively 51 . Groundwater levels in the plain are directly affected by human activities and infrastructure, including industrial and agricultural activities 52 . The climate of the study area is cold and semi-arid. Figure 2 shows www.nature.com/scientificreports/ the precipitation and temperature times series of the region for the period of study. This data was obtained from Tropical Rainfall Measurement Mission (TRMM) satellite precipitation data (https:// disc. gsfc. nasa. gov/ datas ets/ TRMM_ 3B43_7) using the TRMM adjusted merged microwave infrared precipitation rate (in mm/hr) and Root-Mean-Square (RMS) precipitation-error estimates, and Famine Early Warning Systems Network Land Data Assimilation System (FLDAS) datasets (https:// das. gsfc. nasa. gov/ fldas). Times series of precipitation and temperature show average annual precipitation and temperature of ~ 113.8 mm/yr (mostly happening during the rainy season-August to April), and 11.1 °C (with the highest temperature in August and the lowest in February), respectively.
Research methodology. InSAR data. In this study, Sentinel-1 satellite images from October 2014 to January 2021 in a time interval of 12 days were used (Fig. 3). The Interferometric Wideswath (IW) mode was used, which creates images with a 250 km swath at 5 × 20 m spatial resolution. Details of the SAR images that were used for the InSAR time series analysis are presented in Table 1. We processed the LiCSAR generated interferograms over two LiCSAR frames covering Ardabil plain in ascending and descending orbits (LiCSAR frame ID 101A_05193_131313 and 006D_05111_131313) using LiCSBAS package. The interferograms were created with the GAMMA SAR software and filtered using an adaptive phase filter 38,47 . GeoTIFF files of unwrapped phase and coherence are published on the COMET-LiCS web portal on a consistent geographic frame basis in 250 km 2 . Both ascending and descending data were available for the whole or a specific period. The numbers of processed epochs for the ascending and descending frames were 137 and 136, respectively.

LiCSBAS time-series analysis.
LiCSBAS is an open-source package in Python3 to carry out InSAR time series analysis using LiCSAR products (i.e. unwrapped interferograms and coherence) that are freely available on the    36 . LiCSBAS can run all steps selected for analysis with specified parameters at once using a batch script. Figure 4 shows the workflow of LiCSBAS time-series analysis. The first step is to run the script, download LiCSAR products and then convert file format. In the next step, to reduce the impact of unwrapping errors in Ardabil plain, before the subsequent time-series analysis, pixels with a coherence of 0.1 are masked and then the data are unwrapped automatically using SNAPHU 53 . Thereafter, all the unwrapped data and corresponding coherence images are clipped to the Ardabil plain area. The coverage of unwrapped data and coherence threshold are set to < 0.5 and 0.06, respectively. The unwrapped interferograms are resampled and geocoded using the Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM 30 m). The COMET-LiCSAR system processes the InSAR data and generates interferograms that connect each epoch to three or four nearest acquisitions in time, backward and forward. LiCSBAS employs a modified small-baseline NSBAS approach 54,55   www.nature.com/scientificreports/ for the time-series analysis. In the time-series analysis, deficient interferograms are identified and then discarded based on the coherence and coverage of the unwrapped data by checking loop closure 36,56 .
The whole interferometric networks are then inverted for incremental displacements between the acquisition dates, using the least-squares method 36 , to estimate the displacement value of each pixel in the Ardabil plain. Overall, the baseline networks obtained from S1A images for the ascending (A-101) and descending (D-6) tracks include 453 and 365 interferograms, respectively ( Fig. 5a and b). Interferograms that contain noisy data or do not pass a phase loop closure test (indicating severe unwrapping errors) are automatically excluded from the analysis 57,58 . Finally, the maps are high-pass filtered in time and low-pass filtered in space using a Gaussian filter kernel 59 , to separate noise components from the displacement time series (set to 0.50 year-182 days).
In order to compare the LiCSBAS and permanent GPS data, the following equation was used to obtain GPS-LOS measurements in LOS direction 60 : where α is the azimuth of the LOS vector, θ is the incidence angle, and GPS up , GPS e , GPS n are the velocities in vertical directions, the eastwest and northsouth, respectively.

Results
We estimated the surface displacements in the Ardabil plain using LiCSBAS. Then, the accuracy and consistency of the results were analyzed in both LiCSBAS and GMTSAR InSAR tools.
Interpretation of the S1A data by LiCSBAS velocities. The mean velocity maps derived from the time series analysis of S1A images by LiCSBAS are shown in Fig. 6. Several surface displacement signals including subsidence and uplift were detected in the Ardabil plain. Despite having almost all pixels masked in highly vegetated areas due to decorrelation, more than 110 thousand valid pixels (60%), mainly distributed on the Ardabil plain, still remain (Fig. 6). To obtain meaningful and accurate displacement measurements, a suitable reference area, where the displacement signal is almost zero and marked by the black circle in Fig. 6a and b, was considered for the study area. The calculated mean LOS velocities of S1A images from A-101 (measured in the period of November 2014 to September 2020) and D-6 (measured in the period of October 2014 to January 2021) tracks are between −45 mm/yr to + 14 mm/yr and between −53 mm/yr to + 43 mm/yr, respectively. Since the displacement is illustrated in LOS direction, positive and negative velocities represent the motion of the ground towards and away from the satellite, respectively. The mean LOS velocity values show a maximum displacement rate of 45 mm/yr in the LOS direction in an area located in the southeast part of the plain. This area is categorized as a cropland area with excessive groundwater extraction. Results presented in Fig. 6 also show signs of land subsidence in the cropland areas in the northwest region with a maximum displacement rate of 20 mm/yr. Other areas of the plain show no sign of significant subsidence. In Samarin, especially on its side towards Sabalan peak, there is also a large motion of + 43 mm/yr towards the satellite as presented in Fig. 6b. This area is in the mountainous part of the plain and has the lowest alluvial thickness (10-30 m). It is therefore expected to see the lowest rates of land displacement in this area, which is different from our observations in the area. Different factors can cause the observation of such large motion toward the satellite: (1) The measurement of InSAR observations is relative and highly dependent on the choice of the reference area considered for analysis. This reference area is assumed to be a stable region with no displacement. However, in practice, if the reference area is subjected to subsidence, this may cause observation of an uplift in regions with no displacement; (2) The study area is located in a mountainous area, 18 km from Sabalan mountain. In mountainous areas, topography-related tropospheric     www.nature.com/scientificreports/ Figure 8 shows some of the factors that play the most significant roles in land subsidence, including groundwater table, alluvial deposit thickness, type of soil layers and stratigraphy. According to Fig. 8, the southeast part of the plain has the maximum alluvial thickness of 190 m, especially in areas near the cities of Araloyebozorg and Khalilabad. According to studies by Guzy and Malinowska 48 , this area has a clay layer with a thickness of more than 120 m, starting from a depth of about 25 m. The thickness of the sand layer is only 25 m in this area. Also, the ratio of fine-grained to coarse-grained materials in this area is higher compared to the other parts of the plain, increasing potential for land subsidence 63,64 . The groundwater table measurements from the observation wells indicate an average depth of 60-70 m in the southeast, which are deepest in the whole plain. The central part of the plain has the highest illegal and legal well density with significant water withdrawal, and therefore, it is expected to observe higher rates of subsidence in this area compared to other parts of the plain. However, since the groundwater level in this area is high (due to the bowl-shaped nature of the plain), the plain in this area experiences lower rates of subsidence with groundwater withdrawal.

Consistency assessment.
To evaluate consistency of the LiCSBAS package in measuring land subsidence of the Ardabil plain, results of the package were compared with subsidence measurements obtained using the GMTSAR 63 . For this purpose, S1A images of 30 (A-101) and 28 (D-6) tracks, spanning March 2017 to January 2019, were analyzed using the GMTSAR. The interferometric process by GMTSAR based on SBAS-InSAR followed three main steps: First, a preprocessor was used for each satellite data type to convert the native format and orbital information into a generic format. Second, the InSAR processor was used to focus and align stacks of images and convert map topography into phase which is then used in radar coordinates, and form the complex interferogram. The topographic phase component was corrected using the SRTM 30 m DEM. Third, a postprocessor was applied based on GMT, to filter the interferograms and construct interferometric products of   65 . Figure 9 shows the workflow of time-series analysis by GMTSAR. GMTSAR uses an adaptive power spectrum filter and a Gaussian filter for the unwrapping process (set to default filter parameters of alpha = 0.1 and patch size = 32). In general, the numbers of interferograms associated with a perpendicular baseline of smaller than 120 m and a temporal baseline of smaller than 130 days obtained from S1A images of (A-101) and (D-6) tracks are 101 and 97, respectively. Using the least-squares method, the displacement value of each pixel was estimated in the Ardabil plain. To make the results comparable with the LiCSBAS/LiCSAR processing chain, the same multilooking factors of LiCSAR data were employed for GMTSAR processing and the coherence threshold was set to 0.06. Figure 10 compares the LiCSBAS velocities with those measured using GMTSAR for the period of March 2017 to January 2019 for A-101 track. Figure 10a shows the maximum subsidence rates measured using the GMTSAR and Fig. 10b shows the results obtained using LiCSBAS. As observed in this figure, in terms of the maximum subsidence rates, both packages show close values of −53 mm/yr (GMTSAR) and −52 mm/yr (LiCSBAS). For the D-6 track, these values changed to −46 mm/yr (using GMTSAR) and −49 mm/yr (using LiCSBAS), as observed in Fig. 10c and d. Both GMTSAR and LiCSBAS results showed a major land subsidence concentrated in the southeast of the Ardabil plain, near the city of Araloyebozorg and Khalilabad. However, the difference between LiCSBAS and GMTSAR results is the number of masked pixels in areas of low coherence which is significantly larger in GMTSAR.
In Fig. 11, the average coherence maps (March 2017 to January 2019) in Ardabil plain using GMTSAR and COMET-LiCSAR system are presented. Based on the results presented in Fig. 11a and b, the values of average coherence for GMTSAR and COMET-LiCSAR are 0.27 and 0.37, respectively. Both results showed high  Figure 12 shows scatter plots (Fig. 12a) and histograms (Fig. 12b) of average coherence obtained using GMT-SAR and COMET-LiCSAR system. As presented in Fig. 12a for the coherences of 0.15-0.7, the pixels are distributed along the red line, showing consistency between the GMTSAR and COMET-LiCSAR systems. This figure clearly shows a good correlation with an R 2 value of 0.84. Based on the results presented in Fig. 12b, the mean coherence difference between GMTSAR and COMET-LiCSAR is 0.1, with a standard deviation (Std) of 0.01. It is observed that the coherence values produced by GMTSAR are slightly lower than the COMET-LiCSAR system.
Also, for further comparison of differences in coherence values of the two systems (GMTSAR and COMET-LiCSAR), a pair of unwrapped interferogram phases (20170423-20170610) under challenging conditions were analyzed (Fig. 13). As presented in Fig. 13a and b, despite its substantially higher coherence values, the COMET-LiCSAR unwrapped interferogram phases for the scene pair are not different from GMTSAR results for most parts of the Ardabil plain (R 2 = 0.79 and Std = 0.06), and the interferometric phase results have similar patterns. This behavior may be the result of using different approaches to calculate the unwrapped phase using the applied Gaussian filter.
The major difference is in the employed SBAS algorithm time series processing stage. In the implementation of the SBAS algorithm by GMTSAR, whenever there is one NaN value in any of the unwrapped phases, it will skip the whole column before doing the inversion and hence the pixel gets masked out. However, LiCBAS only removes that NaN value from the observation vector and will do the inversion using the remaining valid unwrapped phases. It should also be noted that even in the case of gaps in the small-baseline network (mainly due to decorrelation associated with factors such as vegetation or snow cover), LiCSBAS is able to do the inversion by imposing a linear temporal constraint across pixels with network gap. This explains the larger coverage of the LiCSBAS estimated velocity compared to the one obtained using GMTSAR. Particularly, in areas with relatively lower coherence, LiCSBAS seems to provide better estimate of the land subsidence.
While LiCSBAS is able to accurately detect large-scale land deformation (~ km) of the plain, its accuracy may be considerably affected by the land cover and use type of the area under investigation. Land subsidence measurements of the Ardebil plain obtained using the LiCSBAS and GMTSAR packages were used to evaluate how land use type and cover affect land subsidence estimation using the two packages. For this evaluation, the plain was classified into five different classes based on their landcover types (Fig. 14): Class A is a cropland area that has the highest potential for land subsidence; Class B is the City of Ardebil, classified as urban area with  www.nature.com/scientificreports/ the largest population and located at the central part of the plain; Class C is the airport and its runways, which are located very close to the area with the most subsidence; Classes D and F mainly consist of rangeland and www.nature.com/scientificreports/ dry farming respectively. The displacement rates of these class regions were then estimated using LiCSBAS and compared with those obtained using GMTSAR. According to Table 2, the displacement velocities of LiCSBAS at cropland area for A-101 and D-6 track show subsidence at rates of 52 mm/yr and 47 mm/yr, respectively, which are consistent with GMTSAR. The Std of differences between the derived results in cropland area for both tracks are 0.09 mm/yr and 2.46 mm/yr, respectively. To evaluate consistency between the LiCSBAS and GMTSAR packages, a time series sample of subsidence behavior in Class A region was used. Figure 15 shows the annual displacement time series of an area in Class A region located in the southeast part (yellow triangle in Fig. 14) of the plain measured using LiCSBAS and GMTSAR for Tracks A-101 and D-6 for the period of March 2017 to January 2019. As presented in this figure, for the areas with expected high rates of subsidence, a good agreement is observed between the results of GMTSAR and LiCSBAS. The time series obtained from both packages indicate land subsidence with a relatively high rate of ~ 50 mm/yr. Both analysis packages show a reduction in the subsidence rate in 2019, that can be due to either implementation of some remediation measures including reduction of pumping draft and artificial recharge of aquifers from the land surface or the heavy rainfall in the region in 2019. Figure 16 shows the Tropical Rainfall Measuring Mission (TRMM) and Standardized Palmer Drought Index (SPDI) of the region from 2017 to 2019. The SPDI is a monthly rainfall drought index based on calculation of the monthly cumulative and standardized rainfall anomalies 66 . SPDI is calculated using reference evapotranspiration and precipitation from the 4-km daily GRIDMET (Gridded Surface Meteorological) dataset and a static soil water holding capacity layer https:// devel opers. google. com/ earth-engine/ datas ets/ catal og/ GRIDM ET_ DROUG HT). Dataset presented in Fig. 16 clearly show severe drought and reduced rainfall in 2018 with TRMM and SPDI of In Class B region, as also presented in Table 2, differences in the displacement velocity and Std for A-101 and D-6 tracks are still small. To assess the LiCSBAS and GMTSAR results in Class B region, permanent GPS measurements were used (Fig. 17). The National Cartographic Center (NCC) of Iran (https:// ncc. gov. ir), as a part of the Iranian Permanent Network for Geodynamic (IPGN) program, has established a permanent GPS station (ARDB) near Ardabil. Figure 14 shows the location of the ARDB GPS station by a red triangle. The green arrow indicates the velocity vector with an accuracy of 0.06 mm. Established in 2014, this station has been continuously measuring land deformation since then. Figure 17a shows a behavioral diagram of GPSUp (in vertical direction) measurements that indicates a maximum displacement rate of about 15 mm/yr from 2017 to 2019. The InSAR time-series for the location of the ARDB GPS station along with continuous GPSLOS (in LOS direction) measurements are illustrated in Fig. 17b. Results indicated a very similar subsidence behavior for the analyzed time series, with a root mean square error (RMSE) of around 6 mm/yr. Also, a very good agreement was observed between the InSAR and GPSLOS time-series fluctuations, in terms of both range and data as observed in Fig. 17b.
As presented in Table 2, the LiCSBAS displacement of Class C region for the A-101 track shows subsidence at a rate of 23 mm/yr which is in an acceptable range when compared to the 19 mm/yr subsidence rate measured using GMTSAR. In Class D and F regions, which are the areas with the lowest subsidence rates, the difference between the displacement rates estimated using the two packages were more pronounced. In these regions, displacement rate measurements estimated using LiCSBAS and GMTSAR differed about 6-7 mm/yr for the A-101 and D-6 tracks, respectively.

Discussion
In order to evaluate the subsidence maps obtained from LiCSBAS and their compatibility with real data, it is necessary to adapt it to ground observations. In this study, relationship between land subsidence and groundwater changes, and a geotechnical-based numerical model were used to validate the results of land subsidence obtained using Envisat and S1A images by LiCSBAS.  (Fig. 18b). As a result of these changes in GWL along with low rates of precipitation, great rates of subsidence were recorded that particular area of the plain. According to the results of GWL change, Envisat data (2004-2010) and S1A data (2014-2020) InSAR measurements, every meter drop in GWL has resulted in about 13 mm of annual land subsidence during the years of 2003-2020. Similar observations were reported by 39,64 .
Geotechnical numerical modeling and comparison with SBAS-InSAR results. The land subsidence measurements obtained using SBAS-InSAR results were compared with the results of a rigorous geotechnical numerical analysis of land subsidence performed using FLAC2D in Ardabil plain between the years of 2006 and 2020. FLAC2D (Fast Lagrangian Analysis of Continua in 2 Dimensions) is a numerical modeling software that implements an explicit finite volume formulation (analogous to finite-difference for 2D geometries) to capture the complex behavior of geotechnical systems and soil-structure interaction. For numerical modeling, soil layers were divided into saturated and unsaturated layers based on the location of the groundwater level and then appropriate constitutive models were assigned to each layer to properly capture the stress-deformation behavior of the layers. Terzaghi's 67 one-dimensional consolidation theory was used for the saturated layer and the constitutive model proposed by 67 was used for the unsaturated level. More details about 67 constitutive model were provided in Appendix A 68 . The input parameters for the constitutive models and methods that were implemented for their measurements are presented in Table 3. As presented in Table 3, the input parameters include a series of hydraulic and mechanical properties that were mainly determined by performing a series of simple lab experiments on soil samples taken from the plain. In this www.nature.com/scientificreports/ study, parameters were calibrated based on properties of soil taken from a site near Khalilabad. This area has the lowest GWL and is expected to experience the highest rate of subsidence over the years. The location of the site is presented in Fig. 19. A rotatory wash-boring method was adopted to drill a 105 m borehole in the field. During boring, a SPT (standard penetration test) split sampler was used to take soil samples at 5 m depth intervals. The samples were then packed in rigid boxes and carefully transferred to the advanced soil mechanics laboratory, where index and strength tests were performed on the collected soil samples to measure input parameters for the numerical model. Five Decagon dielectric sensors were also inserted into the drilled borehole to capture and monitor profiles of soil temperature and moisture during an extended time period. The sensor measurements were then used to obtain profiles of suction and calibrate predictions of the numerical model following a procedure described by 69 . Based on the groundwater level measurements presented in Fig. 19, a groundwater level drop of 1 m/yr was assumed for the numerical analysis. Based on field measurements and sensor readings, the groundwater table was assumed to be at the depth of 50 m and a 150 m depth was considered as the depth of the bedrock. Details of soil sampling and model calibration are provided in 15 . Figure 20 compares predictions of the model over the years of 2006-2020 with land subsidence measurements obtained using SBAS-InSAR results. Based on land subsidence measurements, three time periods can be       (Fig. 16b), the plain experienced significant land subsidence. Based on the results obtained using LiCSBAS, During the years of 2018 to 2020, a higher precipitation rate and some controlling measures that were taken by federal agencies decreased the rate of groundwater level drawdown to a value lower than 30 cm/yr. These changes justify lower rates of subsidence measured during the years of 2018-2020. The results of geotechnical numerical modeling showed land subsidence with an average annual rate of 3.8 cm between 2006 and 2020 which was close to what measured using SBAS-InSAR results. Results presented in this figure clearly show a good correlation between model predictions and InSAR measurements with an R 2 value of 0.77 (Fig. 20b). These two parameters are obtained from the consolidation curve in which λ is the normalized consolidation slope in the loading mode and κ is the normal consolidation slope is in the case of loading in saturation condition of ν-ln p' λs & ks λs is the initial wet and dry route slope and ks is the elastic section slope in Wheeler characteristic curve SI*&SD* related to the air inflow (AEV) of the modified SWRC curve of the Wheeler model in two wet and dry modes P' c known as soil pre-consolidation stress, whose value is obtained for consolidation test for different layers for soil saturation, but it should be calculated, for soils that are in the unsaturated state since the hardness or preconsolidation stress of the soil increases by increasing the suction (decreasing the percentage of soil moisture)

k1 & k2
These two parameters are used to model the interaction between mechanical and hydraulic behavior, which moves the LC curve when it comes to SI or SD curves, and vice versa ν0 & γd These two parameters are measurable according to the soil Index characteristics and two parameters of porosity ratio and soil Gs Cν When a time-dependent settlement estimation is required, the coefficient of consolidation should be known. The amount of this parameter can be obtained using the one-dimensional consolidation test (edometer)

Conclusion
In this study, an open-source interferometry time series analysis package, LiCSBAS that integrates with the automated Sentinel-1 InSAR processor was used to measure land subsidence of Ardebil plain, located at northwest, Iran, from October 2014 to January 2021. The results of this study show that LiCSBAS well performed in generating the Sentinel-1 time series over a full six years in the order of nearly 137 processed epochs for the A-101 and D-6 tracks. Using the LiCSBAS analysis package, a mean LOS velocity of 45 mm/yr was estimated for the Ardabil plain. Results also indicated subsidence rates as high as 50 mm/yr in the southeastern part of the plain. From 2018, due to higher rates of precipitation and some prevention measures, lower rates of subsidence were recorded. For consistency assessment of the land subsidence measurements using LiCSBAS, the annual displacement time series of an area near the city of Araloyebozorg (located at the southeast part of the plain) were obtained from March 2017 to January 2019 and then compared with the results of GMTSAR. LiCSBAS provided acceptable results of land subsidence in areas experiencing high rates of displacement and relatively low coherence. Comparison of the results of subsidence estimated using InSAR measurements and those obtained from GPSLOS measurements and a geotechnical based numerical modeling technique showed that LiCSBAS is able to accurately measure land subsidence, especially in areas that experience high rates of subsidence.

Data availability
All the analyzed groundwater level data, climate change and geotechnical data, and GPS observations collected during this study are included in this paper. The InSAR products that support the findings of this study include the following: and ID 006D_05111_131313 include 365 interferograms from (20141006 to 20210114): https:// gws-access. jasmin. ac. uk/ public/ nceo_ geoha zards/ LiCSAR_ produ cts/6/ 006D_ 05310_ 131313/ inter ferog rams/. LiCSAR products based on a predefined LiCSAR frame )from COMET-