Simulation of present and future land subsidence in the Rafsanjan plain, Iran, due to groundwater overexploitation using numerical modeling and InSAR data analysis

In the Rafsanjan plain, Iran, the excessive use of groundwater for pistachio irrigation since the 1960s has led to a severe water level decline as well as land subsidence. In this study, the advantages of InSAR analyses and groundwater flow modeling are combined to improve the understanding of the subsurface processes causing groundwater-related land subsidence in several areas of the region. For this purpose, a calibration scheme for the numerical groundwater model was developed, which simultaneously accounts for hydraulic aquifer parameters and sediment mechanical properties of land subsidence and thus considers the impact of water release from aquifer compaction. Simulation results of past subsidence are calibrated with satellite-based InSAR data and further compared with leveling measurements. Modeling results show that land subsidence in this area occurs predominantly in areas with fine-grained sediments and is therefore only partly dependent on groundwater level decline. During the modeling period from 1960 to 2020, subsidence rates of up to 21 cm year −1 are simulated. Due to the almost solely inelastic compaction of the aquifer, this has already led to an irreversible aquifer storage capacity loss of 8.8 km 3 . Simulation results of future development scenarios indicate that although further land subsidence cannot be avoided, subsidence rates and the associated aquifer storage capacity loss can be reduced by up to 50 and 36%, respectively, by 2050 through the implementation of improved irrigation management for the pistachio orchards.


Introduction
Due to the predominantly arid to semiarid climate and the scarcity of surface water in many parts of the country, groundwater covers about 60% of Iran's water consumption (Mirzaei et al. 2019).Low-cost access to groundwater and energy, regulated by the state, has led to an expansion of the agricultural sector, which is now responsible for approximately 90% of the extracted groundwater (Maghrebi et al. 2020).Also, Rafsanjan's water supply has always depended on the exploitation of groundwater resources (Jaghdani and Brümmer 2010; Goudarzi et al. 2018).From the ancient past until the middle of the last century, this was ensured by shallow hand-dug wells and a distinct qanat system, which provided access to freshwater and thus allowed people to live in the region (Mirnezami et al. 2020).In the late 1950s, people began drilling wells (Mirnezami et al. 2020), becoming less dependent on topography, climatic changes, and natural groundwater fluctuations (Abbasnejad et al. 2016).With the rapidly increasing demand for water, mainly for irrigation of the expanding pistachio plantations in the region, the number of production wells increased from 1961 onwards to such an extent that the aquifer was classified as endangered by the authorities in 1974.As a result, further exploitation and construction of new wells were prohibited (Jaghdani 2012).However, with the decline in regulatory controls beginning in 1980, the number of new legal, but also illegal, wells increased sharply again (Jaghdani 2012).Due to the extensive groundwater extraction, the water level of the Rafsanjan aquifer has declined almost continuously over the past 60 years (Zera'at-kaar and Gol-kaar 2016), resulting in a severe loss of water resources as well as land subsidence, the lowering of the land surface due to a compaction of the aquifer matrix (Motagh et al. 2017).Yet Rafsanjan is not an isolated case.Some of the fastest-sinking cities and regions due to land subsidence are located in Iran (Herrera-García et al. 2021).Due to the aforementioned substantial groundwater consumption for irrigation, it is no surprise that land subsidence occurs predominantly in regions with high agricultural use (Motagh et al. 2008).
In general, subsidence can have various causes and can be observed in many parts of the world.Besides groundwater level decline, land subsidence can be induced by various other natural processes such as tectonics or the dissolution of minerals.However, fluid mobilization in the subsurface due to groundwater depletion constitutes the main driver for subsidence (Herrera-García et al. 2021).In the United States, for instance, it is responsible for more than 80% of the known land subsidence (Galloway et al. 1999).Other prominent examples of land subsidence can be found in, e.g., Bangkok (Phien-wej et al. 2006), Jakarta (Abidin et al. 2011), Mexico City (Osmanoğlu et al. 2011), Tokyo (Hayashi et al. 2009) or California's San Joaquin Valley (Galloway et al. 1999).Typical consequences are the destruction of civil infrastructure such as roads, railways, buildings, pipelines, gravity-driven canals and aqueducts, earth fissures as well as enhanced risk in coastal and river flooding (Hoffmann et al. 2003;Borchers and Carpenter 2014;Haghshenas Haghighi and Motagh 2021).Moreover, the associated compaction of the sediments leads to a reduction in the storage capacity of aquifer systems (Béjar-Pizarro et al. 2017;Herrera-García et al. 2021).
In the past, land subsidence often remained undetected for a long time due to its usually large-scale extent, making a comparison with nearby fixed reference points difficult.This has changed with satellite-based technologies such as the global navigation satellite system (GNSS) or interferometric synthetic aperture radar (InSAR).These methods allow horizontal and vertical terrain deformations to be measured at the Earth's surface with subcentimeter to submillimeter accuracy (Mousavi et al. 2001;Galloway and Burbey 2011;Motagh et al. 2017).These exploration technologies were also applied to monitor land subsidence in Rafsanjan.Installing a global positioning system (GPS) station network, Mousavi et al. (2001) were able to measure the land subsidence in Rafsanjan over 8 months with a vertical and horizontal accuracy of 5-8 and 3-5 mm, respectively.During their measuring period from August 1998 to April 1999, subsidence rates between 1.5 and 12 cm year −1 occurred.Compared to the decline in groundwater level, the ratio of the observed land subsidence is approximately 10% (Mousavi et al. 2001).In addition, Motagh et al. (2017) analyzed land subsidence in Rafsanjan using five different InSAR datasets.The results show subsidence rates higher than 5 cm year −1 , locally even exceeding 30 cm year −1 .These methods are without any doubt very valuable and provide concrete quantitative measurements, but they can only show what has happened in the past and do not describe any underlying processes.A combination of these satellite measurements with numerical groundwater flow models can be very beneficial in this regard.However, using the InSAR analysis to calibrate or validate hydraulic parameters was only applied in a few other studies, e.g., for the Las Vegas Valley in the United States (Burbey and Zhang 2015), the Firenze-Prato-Pistoia in Italy (Ceccatelli et al. 2021), the Saveh basin in Iran (Jafari et al. 2016), and for lowland regions of Semarang City in Indonesia (Lo et al. 2022).Spatially continuous coverage of data provided by satellite measurements can be used to complement local field measurements.Combined with the integration into numerical models, these methods allow a better understanding of hydrodynamic behavior, predictions of groundwater-induced land subsidence, and thus a better assessment of threats to the aquifer system (Hoffmann et al. 2003;Ezquerro et al. 2017).While some studies on groundwater modeling in Rafsanjan plain do exist, they usually focus on socioeconomic aspects, market regulation, and the impact of declining groundwater levels on pistachio production and farmers (Jaghdani and Brümmer 2010; Rahnama and Zamzam 2011;Parsapour-Moghaddam et al. 2015;Babaei and Ketabchi 2020;Akhavan and Gonçalves 2021;Sharghi and Kerachian 2021;Moghaddasi et al. 2022).There are only a few studies that specifically address the modeling of groundwater-induced land subsidence in this region (Toufigh and Sabet 1995;Sayyaf et al. 2014).Nevertheless, they only analyzed past scenarios utilizing just a few monitoring points for local land subsidence.
The aim of this work is, therefore, to combine the advantages of InSAR measurements and numerical groundwater modeling to develop a tool, which allows for the assessment and prediction of land subsidence for the Rafsanjan aquifer system.For this purpose, a transient groundwater flow model is built, integrally coupled with a subsidence module, and calibrated using time series of observed groundwater drawdown and InSAR data.Subsequently, different future scenarios are simulated to illustrate the consequences of groundwater overexploitation.

Study site
The study area, named after the city of Rafsanjan, is located in southeastern Iran in the province of Kerman and has an area of about 4,236 km 2 (Fig. 1).It is composed of the Rafsanjan-Kaboutarkhan plain in the southern part, the Noogh plain in the north-western part, and the Anar-Koshkoueieh plain in the western part (Mousavi et al. 2001).In 2016, the population of the area was 348,111, of which 203,208 lived in urban areas (SCI 2016).
With an annual rainfall of about 100 mm (Motagh et al. 2017) and a potential evapotranspiration of over 3,000 mm (Mehryar et al. 2015), the climate in this region is described as arid according to the Köppen-Geiger classification (Beck et al. 2018).Characterized by long hot summers and cold winters, the climate provides optimal conditions for the cultivation of pistachios.Additionally, subsidies from the government and low costs for water pumping favor production (Akhavan and Gonçalves 2021), making Rafsanjan one of the largest pistachios exporters in the world (Mehryar et al. 2015).High profitability triggered producers in the area to mainly engage in pistachio production (Mehryar et al. 2015), resulting in almost no land being used for other agricultural purposes (Jaghdani 2012).

Geology
The three plains are part of a graben structure (Khamehchiyan et al. 1994), which has been formed by the last Quaternary deformation of the Yazd block, located northwest of Anar (Aghanabati 2004;Ghorbani 2013).With the beginning of subduction of the Neothetic Ocean Plate beneath the Central Iranian Microcontinent in the Eocene, volcanic activity increased, and volcanic rocks were formed.The associated increased tectonic activity is partly responsible for the formation of the graben structure and resulted in a NW-SE orientation of it (Salehi Nejad et al. 2021).As a result of the tectonic events, synforms were formed in the plains as well as antiforms in the mountains (Motagh et al. 2017).In the northeast, the study area is bounded by the Davaran mountain range (Khamehchiyan et al. 1994).The Mesozoic and Early Paleogene rocks found here are mainly limestones and partly clastic sediments (Motagh et al. 2017).In the southwest, the area is bounded by parts of the Urumieh-Dokhtar magmatic belt, which are mainly composed of Tertiary andesitic volcanic rocks (Khamehchiyan et al. 1994).In addition, the Noogh plain is divided from the Anar-Koshkoueieh plain by another antiform.This saddle structure is composed of highly fractured flysch of the Late Cretaceous (Fig. 1).
The plains bounded by the surrounding mountains are mostly built up by alluvial fans that are characterized by coarse-grained sediments such as gravel and sand, thus forming a productive aquifer system (Mousavi et al. 2001).Towards the center of the plains, the Quaternary sediments build thick fine-grained layers which are present in the form of terrace deposits, dunes, and sandy sheets.Due to the high evaporation rates, sabkha sediments such as salt, gypsum and salty silts can also be found (Mountney 2005;Motagh et al. 2017), whereby the latter lead to a decrease in the hydraulic conductivity of the system (Motagh et al. 2017).The bedrock of the study area is composed of limestone and conglomerate towards the northwest (RWCK 2016), and marl near the city of Rafsanjan (Sayyaf et al. 2014).
Based on this information, along with a digital elevation model (NASA JPL 2013), a geological cross section (Khamehchiyan et al. 1994), and 56 borehole logs (RWCK 2016), a three-dimensional (3D) geometric model of the Rafsanjan aquifer system, consisting of four hydrogeological units, was developed.

Hydrology
Mahmoodzadeh and Ketabchi (2021) have summarized various studies to set up a water balance of the Rafsanjan aquifer system.In most cases, the annual groundwater deficit is quite similar varying from 153 × 10 6 m 3 to 201 × 10 6 m 3 .On the other hand, Motagh et al. (2017) estimated a groundwater deficit of about 300 × 10 6 m 3 year −1 between the years 2004 and 2010 based on observed drawdown and InSAR analyses.Furthermore, Babaei and Ketabchi (2022) used the semidistributed WetSpass-M model for the period 2009-2016 to estimate the water balance components, resulting in an average annual groundwater deficit of 256 × 10 6 m 3 .
In this study, the water balance from the Iran Water Resource Management Company (IWRMC 2015) is used as a reference for estimations of the required model boundary conditions because it further subdivides outflows into groundwater extraction, evapotranspiration, and aquifer discharge.In addition, observed groundwater levels for the period 2003 to 2020 as well as detailed groundwater extraction rates for 2010 (IWRMC 2011(IWRMC , 2020) ) are provided by the IWRMC.The evolution of the average groundwater levels from 1951 to 2015 in the aquifer system, derived from Zayandehroodi (2012) and Zera'at-kaar and Gol-kaar ( 2016) is also considered.
In 2010 there were 1,445 official extraction wells in the Rafsanjan district (Zera'at-kaar and Gol-kaar 2016).The primary use of groundwater is irrigation of pistachio fields, which accounted for about 94% of the groundwater withdrawal.The rest of the extracted water is used for industrial and domestic purposes, accounting for about 1 and 5%, respectively (IWRMC 2011).Groundwater extraction rates per well are only available for 2010 (IWRMC 2011); therefore, individual extraction rates per well are estimated based on the temporal evolution of the agricultural area.For this purpose, the annual cultivated area is derived from a land cover classification using images of the Landsat 5 and Landsat 8 missions covering a period from 2003 to 2020.Available images captured in July or August of each year with the least cloud cover were selected.Similar to Mehryar et al. (2015), the land cover classification is computed based on the normalized difference vegetation index (NDVI), performing an object-based classification.Hence, a segmented image is created from the generated raster, which additionally considers closely related pixels.Therefore, areas between pistachio crops are classified as agricultural areas and thus irrigation zones.Subsequently, the individual extraction rates per well for the period from 2003 to 2020 are calculated by multiplying the 2010 rate per unit area by the extent of agricultural land for each year.For the period 1961-2003, total withdrawal rates from different sampling years are used (IWPRI 2012).Data gaps are filled by linear regression-see Fig. S1 of the electronic supplementary material (ESM).

Theoretical background
The vertical compression and expansion of an aquifer are caused by changes in the effective stress acting on the sediment ∆σ′ zz [ML −1 T −2 ], which can result from changes in the total stress ∆σ zz [ML −1 T −2 ], or changes in the fluid pore pressure ∆p [ML −1 T −2 ] (Eq.1; Terzaghi 1925).
Since geohydraulic processes do not cause a change in σ zz (∆σ zz = 0), any change in σ′ zz is attributed to a change in p (Busch et al. 1993).Thus, ∆σ′ zz can be expressed only in terms of the change in the hydraulic head ∆h, the (1) density of the water ρ w , and the acceleration of gravity g (Eq.2).
where db [L] is the change in thickness of a control volume with initial thickness b [L].Accordingly, db represents the compaction (db < 0) or expansion (db > 0) of the control volume.Combining Eqs. ( 2) and (3) yields Eq. ( 4).
S sk [L −1 ] and S k [-] are the skeletal-specific storage and the skeletal storage coefficient, respectively.Note that in the last equivalence of Eq. ( 4), the compaction or expansion only is related to S k and dh.
Depending on whether σ′ zz is lower or higher than the preconsolidation stress σ′ zz(max) , elastic or inelastic deformation occurs within the aquifer system.If σ′ zz equals or exceeds σ′ zz(max) , the deformation is inelastic and it is elastic otherwise (Hoffmann et al. 2003).To consider this, S sk is subdivided into two separate parameters (Eq.5).
where S ske [L −1 ] and S skv [L −1 ] are the elastic and inelastic skeletal specific storage, respectively (Hoffmann et al. 2003).In general, inelastic compression is usually much higher than the elastic one (Hoffmann et al. 2003), being up to three orders of magnitude greater in the case of aquitards (Pavelko 2004).
In the presence of compressible sediments, a second storage term has to be included in the groundwater flow equation to account for storage captures or releases from the compressible interbeds (Eq.6).
where K [LT −1 ] represents the tensor of hydraulic conductivities, W [T −1 ] is the volumetric flux per unit volume of sink/ sources terms, S s [L −1 ] is the specific storage coefficient and S sk [L −1 ] is a representative skeletal storage coefficient for all interbeds in the system.The SUB package computes the storage changes and corresponding compactions for every cell based on Eqs. ( 6) and ( 4), subdividing the storage term in one from the not compressible media and one from the compressible interbeds (Hoffmann et al. 2003).

Observed land subsidence
Considering five different datasets of SAR images from the Envisat, ALOS, and Sentinel-1A satellites, Motagh et al. ( 2017) observed land subsidence rates higher than 5 cm year −1 within an area of approximately 1,000 km 2 of the Rafsanjan aquifer for the periods 2004 to 2010 and 2015 to 2016 using InSAR analysis.In this study, the results of the ascending as well as descending tracks are used as a reference for the calibration of skeletal-specific storages of the model.Land subsidence in the study area is mainly present in two areas: the Noogh Plain and near Hemmatabad-Agah between Koshkoueieh and Rafsanjan with rates locally exceeding 30 and 20 cm year −1 , respectively.To obtain mean subsidence rates for the entire aquifer system, line-of-sight (LOS) deformation velocities derived from the different SAR datasets were converted to vertical displacement rates assuming that horizontal displacement is negligible (Motagh et al. 2017).The vertical displacement rates, hereafter referred to as subsidence rates, were then averaged using a weighting factor according to the temporal extent of each analysis.Since heavy vegetation can lead to coherence losses in InSAR measurements (Galloway and Hoffmann 2007), 1,250 random points were generated which do not intersect with agricultural areas derived from the previously performed land cover classification (Fig. 2).Subsequently, the averaged observed subsidence rates were extracted for these points.To account for uplifting processes, e.g., induced by tectonics, the median of the vertical displacement in areas where no groundwater is extracted (Fig. S2 of the ESM), and hence no subsidence is expected, is used for an offset correction of the subsidence rates.
Other studies have assessed subsidence in the study area as well.Using GPS, Mousavi et al. (2001) measured land subsidence rates of up to 12 cm year −1 , whereas near the city of Rafsanjan rates reached about 6 cm year −1 for the year 1999.Furthermore, Sayyaf et al. (2014) conducted a subsidence simulation using a two-dimensional (2D) finite element model based on topographic leveling campaigns.While leveling data presented vertical displacement of up to 3 cm year −1 in Rafsanjan for the period from 1998 to 2004, their simulated subsidence rates correspond to 4-5 cm year −1 , exhibiting a continuous decrease in the subsidence rate towards the southeast of the city.Although these data are not directly included in the calibration, they still serve as a comparison with the model.
A country-wide InSAR analysis of Sentinel-1 shows that many aquifers in Iran are undergoing long-term subsidence trends with limited seasonal variations due to the recharge of the aquifers (Haghshenas Haghighi and Motagh 2021).Signs of inelastic deformation such as earth fissures and cracks in several areas, suggest that at least parts of the aquifer deformation are inelastic (Motagh et al. 2017).While several studies observed groundwater decline across the country and land subsidence associated with it, no significant rebound has been reported in the past few decades.Therefore, it can be assumed that the groundwater in most of the current subsidence basins is below the pre-consolidation level and that most of the long-term subsidence observed across the country, including in Rafsanjan, is inelastic.On the other hand, in most subsidence areas short-term variations of deformation with small magnitudes compared to the long-term subsidence rate have been observed.Since the short-term variations might be correlated with the discharge/recharge cycle of the aquifers, they can be interpreted as the elastic deformation of the aquifers.To further elaborate on this, an additional line-of-sight (LOS) deformation analysis on the temporal evolution of the subsidence of the Rafsanjan aquifer system was carried out.For this, Sentinel-1 interferograms provided by LiCSAR (Lazecký et al. 2020) were used and processed by the LiCSBAS InSAR time series analysis package (Morishita et al. 2020).

Modeling tools and calibration strategy
The modeling of groundwater flow and land subsidence is performed with MODFLOW-NWT using the Upstream-Weighting Package (UPW) and the Newton Solver (NWT; Niswonger et al. 2011).Land subsidence is simulated by including the integratively coupled SUB Package (Hoffmann et al. 2003).
For the optimization of the hydraulic and subsidence parameters, automatic calibration is used.PEST (Model-Independent Parameter Estimation; Doherty 2015) in combination with the global optimization tool CMA-ES (Covariance Matrix Adaptation Evolution Strategy; Hansen 2006) is applied to minimize the objective function, which is the sum of weighted squared residuals of observed and simulated variables.
The calibration process is divided into two steps.First, the horizontal hydraulic conductivities (K x , K y ) are calibrated during a steady-state simulation, for which the aquifer was subdivided into three plains to account for different sedimentary depositional conditions (Fig. 1).Such a calibration requires data on static water levels, i.e., hydraulic heads that are not subject to any temporal changes induced by pumping.In the absence of these historical data, static water levels were estimated based on the water level evolution (Zera'at-kaar and Gol-kaar 2016).Prior to this calibration, a regularization process was conducted by linking parameters contributing to a null space to corresponding hydrogeologic units.Thus, the null space is avoided and the number of fitting parameters is decreased.For this process, the SSSTAT utility of PEST was used, which calculates statistics of the solution and null spaces based on the singular value decomposition (SVD) of the inverse problem.
In the second step, the vertical hydraulic conductivities (K z ) and the storage terms of S s , S y , S ske, and S skv are calibrated during a transient simulation of the period 1961-2020.These parameters must be calibrated simultaneously because water is released from the aquifer during the simulation of land subsidence, which represents an additional source term in the water balance and thus exerts an influence on the hydraulic system.Before this calibration, also a regularization approach was conducted, this time linking parameters with a null space variance greater than 0.9.To simultaneously calibrate the different storage parameters, two different data sets were used as references.These include data from 78 observation wells spanning a time period from 2003 to 2020 (IWRMC 2020) as well as InSAR-derived subsidence rates between 2004 and 2016 (Motagh et al. 2017).The values of the subsidence rates, derived from the random points (Fig. 2), show a positive-skewed distribution.To focus on areas with significant subsidence rates and likewise avoid the influence of noise due to faults in marginal areas (Motagh et al. 2017), values below the 75%-quantile are weighted with a factor of 0.1 (Fig. S2 of the ESM).A priori parameter ranges of storage parameters for the calibration are derived from Calderhead et al. (2011), Hölting and Coldewey (2013), Chowdhury et al. (2022), andLi et al. (2022).
Since the subsidence-related storage values specified in the SUB package are the elastic skeletal (S ke ) and the inelastic skeletal (S kv ) storage coefficients, the thickness (b) of the different interlayers is not directly considered.Therefore, a calibration scheme (Fig. 3) to account for the aquifer thickness in each cell by decomposing S k into S sk × b was designed.To this end, within an optimization loop, a Python-based routine that multiplies the skeletonspecific storage coefficients generated by PEST by the thickness of the different layers of each cell was added.Thus, discretized skeleton storage coefficients for each cell and layer could be calculated and, subsequently, generate the subsidence input file (*.sub).
Due to significantly different data ranges and units, deviations of the initial heads of each observation well (head), respective drawdown (ddn), and land subsidence (sub) are calculated using scaling factors according to Clausnitzer and Hopmans (2005).This results in the objective function (Eq.7) for determining the optimal parameter combination.
where Φ is the total error and σ is the standard deviation of the respective dataset.
In the transient model, an annual stress period with one time step each is set until the first observed head is available.Afterward, time steps are increased to four per stress period, resulting in a 3-month temporal resolution. (7) n sub × 2 sub Fig. 3 Transient calibration approach.The colors represent the software that is performing the process: gray is PEST, brown is Python, and blue is MODFLOW

Conceptual model
The hydrogeological system of the Rafsanjan aquifer is composed of the three alluvial plains Rafsanjan-Kaboutarkhan, Anar-Koshkoueieh, and Noogh (Fig. 1).The principal groundwater flow direction is from the southeast towards the northwest and hence follows the topography (Jaghdani and Brümmer 2010; Sharafati et al. 2020).Based on the analysis of the lithological profiles, the aquifer system is subdivided into a succession of four hydrogeological units: clayey sand, sandy gravel, silty sand, and gravel, underlain by a calcareous bedrock, assumed as an impervious boundary (Fig. 4).The deepest layer is considered confined, while the upper ones are convertible.The model domain is discretized in a regular horizontal 500 m × 500 m grid.
The annual groundwater recharge rates, resulting from return flows of (1) irrigation and (2) leaking pipes, strongly depend on the amount of extracted groundwater and are calculated considering return flows of 60 and 30%, respectively (Moghaddasi et al. 2022).These rates are distributed according to the principal agricultural and urban areas derived from the land cover classification.In addition, thirdly, there is a constant lateral inflow originating from rainwater infiltrated into tectonic faults and coarse-grained sediments in the surrounding mountains (Motagh et al. 2017), with inflow rates projected according to the area of the respective catchment (Fig. 5a).Fourthly, recharge due to seepage from surface canals, fifthly, diffuse infiltration of precipitation and, sixthly, focused infiltration from surface flow is assumed to be constant over time.Since there is a very widely branched channel network and no permanent river in the study area (Jaghdani 2012), the latter three recharge rates are equally distributed over the modeling area (Table 1).Land classification results for the period from 2003 to 2020 first show an increase in classified pistachio fields, followed by an almost equivalent decrease.Therefore, estimated extraction rates also initially increase from 708 × 10 6 m 3 in 2003 and reach their maximum of 718 × 10 6 m 3 in 2007.With the following decrease in vegetation, the estimated extraction decreases likewise, reaching its minimum of 667 × 10 6 m 3 at the end of the modeling period in 2020 (Fig. 5b).The total groundwater extraction for the entire modeling period is presented in Fig. S1 of the ESM.
Two constant heads are defined as outflow at the northern boundaries of the aquifer system (Fig. 5a), with starting and ending heads derived from the temporal evolution of the nearest existing observation well.In the case of the steady-state simulation, static water levels are extrapolated based on the average groundwater decline in the period of 1960-2010 (Fig. 5c).
The time series of the lowering of the land surface in Rafsanjan confirm that the subsidence rates in different years do not follow rainfall patterns, which can result in higher groundwater recharge rates, hence leading to a recovery of the groundwater levels.Consequently, if the land deformation is elastic, the increase in water pore pressure would induce an uplift of the land surface (Wilson and Gorelick 1996).Figure 6 presents the results of two exemplary time series for points near Hemmatabad and Noogh (locations provided in Fig. 2), showing primarily long-term subsidence trends with small magnitudes of short-term variations.The time series near Noogh represents a significant linear trend of subsidence of as much as 180 cm from November 2014 to September 2022.The other time series, however, shows as much as 150 cm subsidence, with signs of deceleration throughout the 8 years.Furthermore, the time series of deformation plotted for both points do not show a significant correlation with precipitation-for example, neither higher than average precipitation in 2017 and 2019 nor lower than average precipitation in 2018 and 2021 are followed by acceleration/deceleration in subsidence.Assuming that groundwater extraction rates have been almost uniform in the past eight years, the results confirm that the deformation is mainly inelastic.

Hydraulics
Through the regularization process, 29 parameters turned out to be not very sensitive and consequently tied during the calibration.On the other hand, the inelastic skeletal specific storage coefficient (S skv ) appeared to be particularly relevant; therefore, only 2 of the 12 S skv parameters were tied after regularization, uniting the sandy gravel layer under a single parameter.The complete list of tied parameters can be found in Table S1 of the ESM.All of the calibrated parameters show sedimentcharacteristical values for the respective layers (Table 2).While mixtures of silty sand with clay interbeds are generally classified as sandy silt or silty sand, depending on their hydraulic conductivities, gravelly layers are at the transition between sandy gravel and gravelly sand (Busch and Luckner 1974).Nevertheless, calibrated hydraulic conductivities may be higher than literature values due to the occurrence of tectonic faults in the study area (Altunkaynak and Şen 2011; Motagh et al. 2017).Calibrated inelastic skeletal specific storage for clayey and silty sand layers correspond with values derived in different studies (Hoffmann et al. 2003;Rezaei et al. 2020;Li et al. 2022).In agreement with Sneed (2001), calibrated S skv of the sandy gravel layer is approximately two orders of magnitude lower than in the finer-grained layers.
Overall, the model exhibits a water balance error of up to 0.010% with an average of 0.001%.Groundwater released from aquifer storage peaked at 474 × 10 6 m 3 year −1 in 1989 (Fig. 7).Afterwards, it slightly decreased and stabilized at 423 ± 15 × 10 6 m 3 year −1 until 2020.The cumulative water released from aquifer storage reached 20.5 km 3 by 2020, of which 8.8 km 3 (i.e., 43%) correspond to subsidence-induced water release.Given that the subsidence in the Rafsanjan aquifer occurs predominantly by inelastic compression, it can be assumed that these 8.8 km 3 correspond to a permanent loss of the aquifer's storage capacity.
For the steady-state model's calibration, an NRMSE (normalized by the difference between the highest and lowest observation) of 0.076 is obtained.With a PBIAS of -0.9%, the model shows a minimal underestimation of the initial heads (Fig. 8).Higher head residuals near the transition of the different plains (Fig. S3 of the ESM) might be caused by the abrupt change of the hydraulic conductivities of each zone representing different sedimentary deposition conditions.These changes are more likely to occur in a continuous course under natural conditions.In the transient model, the fitting of the initial heads of each observation well shows an NRMSE of 0.083, and a slight underestimation of the heads with a PBIAS of -0.2% (Fig. S4 of the ESM).
Groundwater drawdowns are one of the important aspects of land subsidence modeling since the behaviors are directly proportional (Eq.4).In the study area, two significant centers of high water level decline with a total drawdown of  up to 35 m are simulated, whose location correlates as expected with the pumping rates, i.e., higher drawdowns in agricultural areas than in urban areas (Fig. 9).On the other hand, the lowest drawdowns are located in the northwestern boundary of the Noogh plain, where observed water levels are often nearly constant.
Since the distribution and percentage of extraction of wells over the entire modeling period were assumed constant, potential unknown changes in well locations due to relocations of pistachio crops to more efficient areas (Mehryar et al. 2016) may also lead to under-or overestimation of pumping rates in different areas.As a result, slight differences between the simulated and the first measured head of each observation well can occur in the model (Fig. S4 of the ESM).However, the drawdowns occur predominantly in the corresponding layers and display a quite good fit in most of the observation wells (Fig. 9), resulting in an overall RMSE of 3.49 m, an NRMSE of 0.139, and a PBIAS of -11.2%.

Subsidence
Interestingly, the simulated land subsidence (Fig. 10a) does not necessarily correlate with the groundwater drawdown, as low or no land subsidence rates are simulated in several areas where large depression cones are present.Therefore, land subsidence in Rafsanjan plain is not solely caused by groundwater depletion.Instead, the occurrence of the highest land subsidence tends to be in areas with thicker layers of finer-grained sediments due to the possible nonreversible rearrangement of grains resulting in an inelastic compaction (Wilson and Gorelick 1996;Li et al. 2023).
The results of the subsidence simulation exhibit a RMSE of 2.5 cm year −1 , an NRMSE of 6.7%, and a PBIAS of 1.4%.This slight overestimation might be due to the simulation of subsidence rates in large areas where nearly no subsidence has been observed (<75%-Quantile in Fig. 10b).
The highest subsidence rates are evident within the Noogh plain, where also thicker layers of fine-grained sediments like clayey sand and silty sand are present, even though the observed water level is relatively small (Fig. 9).These areas show the highest skeletal storage coefficients (S k ; Fig. S3 of the ESM).Subsidence between the pistachio orchards in this region is slightly overestimated, while higher subsidence rates in the crop areas are, in general, slightly underestimated.Therefore, the model reproduces an average aggregated subsidence within the Noogh plain.The second center of land subsidence between Koshkoueieh and Rafsanjan is shifted slightly towards the northwest in comparison to the depression cones due to the distribution of the thickness of fine-grained layers.In the areas outside the subsidence centers, predominantly low S k values correlating with low subsidence rates are observed.While southwest of Rafsanjan, subsidence rates are slightly underestimated in the model, simulated land subsidence rates are higher than observed in the north of Koshkoueieh.This overestimation might be a result of the interpolation of the geometry that is based on only one borehole in this area.Additional misleading interpolation is suspected in the southeastern part of the study area since high subsidence rates are simulated here, although almost no subsidence has been observed.
Nearly no subsidence is simulated in the marginal areas of the aquifer, where predominately coarse-grained sediments are present.InSAR measurements in these regions present very low land subsidence, and sometimes even uplift.However, variations in the observed data can be caused by changes in subsidence rates during different measuring intervals and due to slight differences in the angles of incidence and heading between different sensors (Motagh et al. 2017).Ranges of the respective data sets are shown as error bars in Fig. 10b.
In comparison to Mousavi et al. (2001), slightly higher subsidence rates were simulated (up to 20.1 cm year −1 ) than the ones measured by GPS surveys over the entire study area between 1998 and 1999 (up to 12 cm year −1 ).However, the highest subsidence rates occur within the Noogh plain, where the density of GPS stations is very low; hence, higher subsidence rates could probably have been observed if the observation density within the subsidence centers had been increased.Near the city of Rafsanjan, subsidence rates of up to 5.3 cm year −1 are simulated, slightly below the 6 cm year −1 derived from GPS measurements.
In the southeastern part, along the Rafsanjan-Kaboutarkhan road (Fig. 2), the model predicts nearly the same  (2014).Near the city of Rafsanjan, subsidence rates of 3.5 and 4.5 cm year −1 are measured and simulated, respectively, while the model simulates subsidence rates of 4.3 cm year −1 .Table 3 shows a comparison between modeled subsidence and the ones observed in previous studies.

Scenario analysis
Three different future scenarios from 2020 to 2050 are created to assess the impacts of pumping on subsidence in the study area.For all scenarios, anthropogenic groundwater recharges remain dependent on groundwater abstracted, while natural groundwater recharges stay constant.During the first scenario, it is assumed that there will be no change in groundwater extraction.The second and third scenario consider a decrease in groundwater extraction due to improvements in irrigation management.Currently, pistachio fields in Rafsanjan are managed almost exclusively by surface irrigation methods such as flood and furrow irrigation (Sedaghat 2008;Sedaghati et al. 2012).By using more modern methods, like subsurface or micro irrigation, and improved irrigation scheduling, water used for pistachio orchards can be significantly reduced (Goldhamer 2005;Iniesta et al. 2008).Comparing these irrigation systems to traditional methods, Sedaghati et al. (2012) observed no significant changes applying only 60% of the original irrigation demand, while dropping it to 40% resulted in a detrimental effect on pistachio yields.Therefore, by introducing and combining different methods, possible savings of 30% and Unchanged groundwater extraction from 2020 to 2050 leads to total subsidence of up to 5.6 m within the Noogh plain and 4.8 m between Koshkoueieh and Rafsanjan.Accordingly, average subsidence rates of up to 18.7 cm year −1 will be reached (Fig. 11a).With a 30% decrease in groundwater extraction, land subsidence of 3.6 m is expected to occur within the Noogh plain as well as around the region of Hemmatabad-Agah (Fig. 11b).Moreover, reducing groundwater extraction by 50% can effectively minimize further subsidence within the Noogh plain to a maximum of 2.8 m from 2020 to 2050.Interestingly, in this scenario, the highest land subsidence of up to 3.1 m (Fig. 11c) takes place between Rafsanjan and Koshkoueieh.This shift of the subsidence hotspot is attributed to the different layers with different S skv being affected by changes in extraction rates and corresponding drawdown over time.However, within the city of Rafsanjan, the total subsidence during the scenario simulation period would be less than 1.0 m if the groundwater extraction is reduced by 50% (Fig. 11c).Loss of aquifer storage capacity due to land subsidence cannot be avoided in any scenario.In the case of constant groundwater extraction, an additional 5.9 km 3 of storage capacity would be lost until 2050.Reductions in groundwater extraction by 30 and 50% resulted in an irreversible loss of aquifer storage capacity of 4.5 and 3.8 km 3 , respectively.

Conclusions
Using a combination of satellite-based methods and a hydrogeological model, land subsidence due to groundwater overexploitation was simulated in the Rafsanjan aquifer system.Compared to InSAR analysis, the simulated values exhibit an RMSE of 2.5 cm year −1 (NRMSE of 6.7%).While satellite-based methods provide a good general overview of the S k of the entire sedimentary succession of an aquifer system, groundwater models allow accounting for partial vertical differences.In this study, a calibration scheme, which enables considering the cell-dependent aquifer thickness by decomposing S k into S sk × b was designed and, additionally, accounts for the interaction between land subsidence and hydraulics.The results show that the release of water from land subsidence has a significant impact on potentiometric heads and the water budget, with 43% of the total withdrawal from the aquifer storage originating from subsidenceinduced groundwater release.Moreover, it has been shown that for the Rafsanjan aquifer system, material properties of sedimentary layers have a greater impact on land subsidence than the magnitude of groundwater level decline.Finally, the model was applied in a forward simulation of different development scenarios.It can thus be predicted that constant groundwater extraction until 2050 will lead to further land subsidence of up to 5.6 m in some areas, causing damage to technical infrastructures such as roads, houses, and water supply systems.Furthermore, groundwater extraction leads to irreversible changes in the aquifer system, e.g., it already resulted in a loss of aquifer storage capacity of 8.8 km 3 and is projected to become 14.8 km 3 in 2050, assuming a continuous abstraction of groundwater.To minimize these damages, water consumption for pistachios should be significantly reduced in the future by improving irrigation management.With a decrease of 50% in groundwater withdrawal used for pistachio irrigation, land subsidence can maintain fewer than 1 m in the city of Rafsanjan, and the aquifer storage capacity loss can be reduced by up to 36% until 2050 compared to the status quo scenario, assuming unchanged abstraction rates.In addition to irrigation management adjustments, it would also be beneficial to implement ground-based geodetic monitoring of land subsidence.Such data could complement model predictions and help to further improve the robustness of the model.DEAL.

Fig. 2
Fig. 2 Subsidence rates derived from Motagh et al. (2017), leveling stations (green squares) derived from Sayyaf et al. (2014), locations of line-ofsight (LOS) deformation time series (yellow stars), and randomly generated observation points for the extraction of subsidence rates derived from InSAR analysis

Fig. 4
Fig. 4 Geological model of the Rafsanjan aquifer system

Fig. 5 a
Fig. 5 a Conceptual model presenting observation wells (IWRMC 2020), lateral inflow contributing watersheds, and extraction rates of 2010 (IWRMC 2011); b Pumping rates derived from land cover classification; c Average drawdown in the aquifer

Fig. 8
Fig. 8 Calibration results for groundwater heads under steady-state conditions

Fig. 9
Fig. 9 Drawdowns relative to steady-state condition (1961-2020).The red y-axes correspond to the observed and simulated water levels in m a.s.l.

Fig. 10 a
Fig. 10 a Observed and simulated land subsidence; b Comparison of residuals with the weighted RMSE of the total dataset

Table 2
Calibrated parametersa Refer to Eq. (6) b Simulated as a confined layer, therefore S y was not considered