Has the threat of groundwater subsidence in Iran been revealed by examining the vast extent of irreversible groundwater depletion by analyzing the land surface across the country of Iran? Why has Iran resorted to groundwater depletion instead of managing watersheds, building dams on rivers, and creating forests and artificial lakes, causing groundwater subsidence and destroying its environment? Is there watershed management in Iran?

Abbas Kashani

University of Mohaghegh Ardabili

{Jessica A. Payne, A. R. Watson, Y. Maghsoudi, S. K. Abmeyer, R. Rigby, M. Lazki, M. Thomas, and J. R. Elliott}

Abstract Ongoing depletion of Iran's groundwater, driven by human extraction, has contributed to 106 incidences of basin‐scale, land‐surface subsidence covering 31,400 km2 (>10 mm/yr, 1.9%) of the country. We use Sentinel‐1 Interferometric Synthetic Aperture Radar time series to map and analyze, for the first time, surface velocities within these subsiding regions in vertical and east‐west directions. We find maximum subsidence rates in the vertical direction reach ∼340 mm/yr in Rafsanjan, with 77% of subsidence faster than 10 mm/yr correlating with agriculture. We assess the risk posed by differential subsidence to residential populations, estimating that ∼650,200 people in Iran are exposed to medium or higher subsidence induced risk caused by steep differential subsidence gradients. We further demonstrate the use of these vertical and east‐west velocity gradients in aiding identification of structural and geological controls on subsidence patterns, some of which are not evident on existing fault maps. We use Independent Component Analysis nationwide to separate subsidence deformation sources and demonstrate that most of Iran's rapid subsidence, and thus aquifer compaction, is irreversible, with inelastic deformation contributing at minimum 60% of the observed deformation magnitude. This proportion of deformation which is irreversible varies within and between subsidence regions. During a recent, severe regional drought (2020–2023), we demonstrate the control of precipitation on the elastic, recoverable subsidence deformation magnitude, with the elastic to inelastic deformation ratio falling from 41% to 44% pre‐drought to 31%–36%.

1. Introduction

Billions of people and nearly half of all irrigated agriculture rely on groundwater as a primary water source (Famiglietti & Ferguson, 2021). Groundwater provides ∼60% of Iran's total water supply, with 90% of this groundwater withdrawal supplying agriculture (Noori et al., 2021). Iran has for at least two and a half millennia used systems of underground qanats (tunnel systems) for groundwater transfer to the surface (English, 1968; Noori et al., 2021). In more recent decades, technological advances and modernization of agriculture caused qanats to dry out due to deepening groundwater levels (Noori et al., 2021). Such advances followed rapid, nationwide urbanization beginning in the 1950–60s (Mashayekhi, 2016), which transformed agricultural land into urban and industrial land in some regions (Maghrebi et al., 2021), and expanded agricultural land elsewhere in the country. During this expansion period, higher than average precipitation in the early 1990s (Madani et al., 2016) led to the construction of several large dams (Alizadeh & Keshavarz, 2005), which contributed to a two‐thirds

M. Thomas, J. R. Elliott

Plain Language Summary Stores of underground fresh water are shrinking in Iran as people are removing this water to grow crops. When these stores shrink, this can cause the land surface to move down. We use satellite images to find 106 places where the surface in Iran is moving down quickly due to shrinking water stores. More than 31,400 km2 of Iran (an area the size of Belgium) is moving down faster than 10 mm/yr. Threequarters of this land is used for farming or buildings. We can use this data to tell us where cracks or faults (where earthquakes can happen) might be hidden under the land surface. Karaj has the highest chance of these cracks causing damage to buildings—23,000 people live where there is high chance of damage due to downwards land movement around these cracks. We find most of this land movement due to water store shrinking cannot be reversed. Recently (2020–2023) there was less rain in Iran, which meant even more downwards land movement was not able to be reversed. Finally, sometimes we see that the surface can actually move upwards around where the land moves down.

increase in irrigated farmland area. During Iran‐wide extreme droughts of the late 1990s and 2008, the use of technologies such as deep well drilling, and surface water damming and diversion increased to maintain groundwater supplies and food security. This increase contributed to the number of groundwater extraction points increasing ∼85% 2002–2015 (Noori et al., 2021). Despite these strategies, annual groundwater withdrawal still decreased 18% 2002–2015 to 61.3 km3/yr, in part due to depletion or salinization of fresh groundwater resources (Noori et al., 2021). This declining recharge alongside non‐renewable groundwater extraction has led to 77% of Iran's land area experiencing extreme overdrafts (when groundwater pumping rates exceed recharge rates; Ashraf et al., 2021; Jasechko et al., 2024). This groundwater stress is not equal across Iran, with river basins around Tehran, for example, experiencing extreme water scarcity whilst south‐east Iran experiences mild water stress (Müller Schmied et al., 2021). An important consequence of groundwater depletion, and the focus of this study, is land surface subsidence: the lowering of the Earth's surface that mainly occurs due to solid or fluid mobilization underground (Herrera‐García et al., 2021). This process is driven by reduced pore water pressure (Section S2 in Supporting Information S1, Terzaghi, 1925; Leake, 1990; Galloway et al., 1998).

Due to deepening groundwater levels and extreme overdrafts (Ashraf et al., 2021; Jasechko et al., 2024; Noori et al., 2021), Iran experiences some of the fastest land surface subsidence rates in the world, with peak subsidence rates exceeding ∼200 mm/yr in some regions (Haghshenas Haghighi & Motagh, 2024; Mahmoudpour et al., 2013; Mirzadeh et al., 2021; Motagh et al., 2008). Land subsidence itself presents a hazard to buildings and infrastructure, particularly those located on steep gradients of surface subsidence velocity, or spanning boundaries between differential subsidence rates (Cigna & Tapete, 2021b; Fernández‐Torres et al., 2022; Haghshenas Haghighi & Motagh, 2024; Solano Rojas et al., 2017). The impacts of differential subsidence are observed across several regions of Iran including around Tehran (Mahmoudpour et al., 2013) and Rafsanjan (Motagh et al., 2017). In regions of rapid, urban subsidence elsewhere in the world, studies have used differential settlement derived from Interferometric Synthetic Aperture Radar (InSAR) to map hazard due to subsidence induced fissures and bedrock faults—particularly in Mexico—using angular distortion (Cabral‐Cano et al., 2008, 2015; Castellazzi et al., 2017; Cigna & Tapete, 2021b; Fernández‐Torres et al., 2022). Knowing the location, extent, and differential rates of subsidence is necessary to mitigate against the impacts of land subsidence through improved and informed groundwater management.

Traditionally, land subsidence rates and extents are measured using levelling, extensometer, and Global Navigation Satellite Systems (GNSS) data. Data measured using these techniques are often geographically sparse point data (Chaussard et al., 2013; Herring, 1999), with all three techniques often requiring site visits which can be expensive and infrequent. For GNSS stations that record motion in three components, the vertical component can be used to analyze subsidence rates over time (Baldi et al., 2009; Susilo et al., 2023). However, GNSS vertical velocity uncertainties are often systematically 2.5–3 times greater than horizontal velocity uncertainties, regardless of velocity magnitudes (Sokhadze et al., 2018), despite GNSS velocities being in a global reference field. Moreover, only a few of these three‐component stations record non‐tectonic subsidence in Iran, as many Iranian stations were installed by the Iranian Permanent Network for Geodynamics (IPGN) with deep rooted concrete pillars on stable ground to study tectonic behavior (Khorrami et al., 2019).

Alternatively, InSAR time series analysis provides a complementary approach to GNSS, allowing land subsidence velocity mapping over extensive geographical areas with high precision and improved spatial resolution (Chaussard et al., 2013). Additionally, InSAR vertical velocity uncertainties—as estimated using velocity decomposition methods such as in this study—are systematically half the magnitude of east‐west InSAR velocity uncertainties (Figure S6 in Supporting Information S1), although they are in a local reference field. Due to these advantages, many studies over the last 20 years have mapped and analyzed subsidence in Iran using InSAR. Haghshenas Haghighi and Motagh (2024) conducted the most widespread assessment of subsidence in Iran, using descending Sentinel‐1 InSAR data to map line‐of‐site (LoS) time series velocities across Iran. These authors find that over 40,000 km2 of the country is subsiding faster than 10 mm/yr. Other studies calculate subsidence velocity fields for one or a few regions in Iran using other C‐band or L‐band SAR satellite data (Anderssohn et al., 2008; Foroughnia et al., 2019; Mahmoudpour et al., 2013; Mirzadeh et al., 2021; Motagh et al., 2008).

Subsidence within aquifers results from primary consolidation due to subsurface fluid depletion and secondary consolidation due to geological timescale overburden pressure (Liu et al., 2020). Due to the geologically short observation period of Sentinel‐1 and other satellites, only primary consolidation is likely to impact satelliteobserved subsidence. Primary consolidation consists of both permanent aquifer compaction, which mainly.

📷

Figure 1. Plains of Iran from Maghrebi et al. (2021); Noori et al. (2021) as of 2018. Free plains (215) are where permits are issued for drilling new wells. Critical plains (31) are projected to become prohibited plains in the future. Critical prohibited and prohibited plains are areas where new drilling is banned except for potable water (127 and 236 respectively). In grey are the 106 regions where subsidence is observed with Interferometric Synthetic Aperture Radar (this study) to subside faster than 10 mm/yr over an area larger than 4 km2. Some correlation between the two data sets is observed. The 10 most populated cities in Iran representing 22.5 million people (∼27% of the population; SCI, 2016) are labelled.

occurs in aquitards (inelastic deformation), as well as recoverable compaction which could be reversed through improved groundwater management (elastic deformation; Chaussard & Farr, 2019; Mirzadeh et al., 2021, 2023). Independent Component Analysis (ICA) is a method of blind signal separation which has previously been used to extract these elastic and inelastic deformation components from long‐term, land subsidence InSAR time series (Chaussard et al., 2017). These isolated components have been used to demonstrate the reduced magnitude of elastic recovery of aquifer compaction during drought periods (Chaussard & Farr, 2019).

Here we expand on this literature to use eight years (2014–2022) of ascending and descending Sentinel‐1 radar acquisitions combined with horizontal GNSS station velocities to derive surface velocity fields in two component directions (vertical and east‐west) for 106 rapidly subsiding regions across Iran (Figure 1) at ∼100 m spatial resolution. We estimate and map subsidence induced risk for highly populated cities, as well as demonstrate the use of spatial velocity gradients angular distortion β and horizontal strain ϵ in providing evidence for geological structures such as intra‐basin faults which are not apparent on existing maps. We use blind‐source separation technique ICA to automatically and manually isolate inelastic and elastic components of deformation within subsiding regions, respectively. Finally, we find that the ratio of elastic to inelastic deformation in the Western Tehran Plain reduces during a regional drought, demonstrating the reduced ability of aquifer storage to recover elastically during prolonged dry periods. Our findings are of relevance to groundwater management approaches in Iran, such as informing drilling plain status maps (Figure 1).

2. Methodology

2.1. Region Selection

Subsiding regions are identified based on several criteria. Firstly, candidate regions must be subsiding faster than 10 mm/yr (on average over the past decade) over a geographic area exceeding 4 km2, as measured in our data. We chose this velocity threshold to discount regions apparently subsiding due to accumulated phase bias (e.g., De Zan et al., 2015) or other noise signals. We test the magnitude of phase bias in our subsidence velocities (Maghsoudi et al., 2022) and find that phase bias contributes < ∼ 5% uncertainty (Figure S1 in Supporting Information S1).

Candidate subsiding regions must additionally be used by humans, be it for agriculture or settlements, as we assume groundwater is pumped and used in the same or nearby vicinity. This criteria discounts regions that might accumulate a fast, non‐aquifer depletion related vertical signal. For example, we exclude subsidence near coastal cities Babol and Rasht in northern Iran due to uncertainty of the drivers of land subsidence in these locations, as well as noise contributions to velocities and low coherence due to these region's marsh land cover (Figure S2, Section S1.1.1 in Supporting Information S1).

Finally, we exclude potential candidate regions where the identified subsidence signal corresponds to an earthquake exceeding MW 4 in the observation period (2014–2022).

We name subsiding regions using either names from the literature, or an English translation of the Farsi plain name (one of ∼600 units in Iran defined by a common watershed) in which the subsiding region is located (e.g., Noori et al., 2021). If the region spans two plains then both plain names are used. If two subsiding regions are located in the same plain then they are numbered, for example, Bardsir 1 and Bardsir 2 (Table S1; Figure S5 in Supporting Information S1).

2.2. Short Baseline Interferometric Synthetic Aperture Radar (InSAR) Time Series Analysis

Time series analysis with Looking into the Continents from Space (LiCS) interferogams (Lazecký et al., 2020) and LiCSBAS software (Figure 2, Figure S3 in Supporting Information S1; Morishita et al., 2020) is used to estimate pixel‐wise line‐of‐sight (LoS) cumulative deformation and surface velocities. Generally, we do not downsample any ∼100 m LiCS interferograms during time series analysis so as to obtain more accurate locations of subsidence margins and interior structures. For consistency, we do downsample, however, interferograms in four frames which have a finer pixel width (∼50 m) on the LiCS portal (Lazecký et al., 2020; Table S2 and Figure S3 in Supporting Information S1).

We apply atmospheric noise correction to each epoch using the Generic Atmospheric Correction Online Service for InSAR (GACOS; Yu et al., 2017; Yu, Li, & Penna, 2018; Yu, Li, Penna, & Crippa, 2018) as per methodology of Watson et al. (2024).

To improve the quality of input interferograms for time series analysis, we choose to mask pixels with coherence less than 0.1. Masked interferograms are then clipped to a rectangular extent slightly larger than each subsiding region to ensure stable ground is captured in each image as a zero reference level for deformation. Masked, clipped interferograms with low average coherence (less than 0.05) or with unwrapped data covering less than 60% of the LiCSAR frame are then discarded to improve the quality of the LoS inversion.

Before inversion, we check for unwrapping errors using phase misclosure (Biggs et al., 2007) and remove unwrapped pixels from the time series if any loop that pixel is part of has Φ123 exceeding |π|. This somewhat destructive approach to mitigating loop closure errors is necessary because most subsiding regions in Iran are covered by farmland (Figures 5e and 5f) and are therefore likely to be vegetated. These vegetated regions are both more likely to become incoherent and also to induce unwrapping errors. Reducing loop closure errors in this aggressive manner is one approach we use to mitigate against incoherence induced unwrapping errors. The pixel in each frame with the lowest root‐mean‐square of all loop phase closures which that pixel is involved in is selected as the reference point (Morishita et al., 2020).

After these quality checks, LoS displacements are estimated and referenced following interferogram inversion as per Morishita et al. (2020). Velocity uncertainties are estimated following Watson et al. (2024). Finally, LoS velocities are masked using noise index thresholds defined in Table S3 in Supporting Information S1. Further explanation of each index can be found in Morishita et al. (2020). For frames with interferograms with persistent unwrapping errors, we apply a reunwrapping approach detailed in Lazecký et al. (2024).

2.3. Line‐Of‐Sight Velocity Decomposition

To calculate subsidence magnitudes and rates in the vertical direction due to groundwater depletion, and to aid the location of faults and fissures using horizontal surface motion, we decompose LoS velocities into vertical and horizontal (east‐west) components (Figures 2 and 3; Ou et al., 2022; Watson et al., 2022) as per methodology of Watson et al. (2024). Where subsiding regions exceed the boundaries of a LiCS frame in the along‐track direction, we first merge along‐track adjacent frames' LoS velocities. We correct for plate motion in the LoS direction.

📷

Figure 2. Flowchart visualizing the data‐related methodology used in this study. Parallelograms represent data with grey parallelograms indicating an output data set; no shading represents an input data set. Rectangle boxes indicate a stage of processing conducted in this study. Question marks indicate that not all LoS frames' components include an obvious elastic signal of deformation when compared to precipitation time series. COMET = Centre for Observing and Modelling Earthquakes, Volcanoes and Tectonics; LiCS = Looking into the Continents from Space; UNW = unwrapped interferograms; COH = coherence; A = ascending; D = descending; LoS = line‐of‐sight; ICA = Independent Component Analysis; VU = vertical velocities; σU = vertical velocity uncertainties; VE = east‐west velocities; σE = east‐west velocity uncertainties; β = angular distortion; ϵ = horizontal strain; dU = absolute maximum vertical displacement gradient in space; dE = total horizontal displacement in the east‐west direction; l = pixel width; ppt = precipitation; t‐s = time series.

following the methodology of Stephenson et al. (2022). We then reference LoS velocities for each merged track into a Eurasia‐fixed reference frame defined by an interpolated 2‐D, horizontal GNSS field using a planar, first order polynomial fit whilst masking LoS velocities ±10 mm/yr. We do not use the filtering‐based approach to GNSSInSAR referencing used by Watson et al. (2024), although we do use the same interpolated GNSS data set (based on Khorrami et al., 2019). There is therefore no vertical GNSS used during referencing. During LoS decomposition we decompose overlapping subsiding regions together to reduce uncertainty. For example, regions Isfahan ‐ Barkhar, Lenjanat, and Kouhpaye ‐ Sagzi in central Iran are decomposed from LoS to component velocities in the same inversion. Where these regions overlap, the resulting uncertainty in both components is reduced (Figure S6 in Supporting Information S1).

For frames 028A_05385_191813 and 035D_05397_131013 (Tehran and the surrounding region, Figure S4 in Supporting Information S1) we repeat the above time series procedure at ∼50 m pixel width.

2.4. Comparison to Potential Subsidence Hazard

We compare the locations of measured subsidence in Iran (vertical velocities ≥10 mm/yr over an area larger than 4 km2) to those predicted by Herrera‐García et al. (2021). This study conducted a statistical analysis of lithology, land‐surface slope, land cover, and climate classes to estimate global subsidence susceptibility. By aggregating this susceptibility with the probability of groundwater depletion they estimated potential global subsidence at 30 arc‐second spatial resolution.

📷

Figure 3. Example input data and velocity decomposition results for subsidence region Semnan (Figure S5 in Supporting

Information S1) in north‐central Iran. (a) ).

In previous literature, both population density (Cabral‐Cano et al., 2015; Fernández‐Torres et al., 2020) and building density (Cigna & Tapete, 2021a; Fernández‐Torres et al., 2022) have been used to estimate exposure of buildings to subsidence induced hazard (Cigna & Tapete, 2021a). In Iran, there is a lack of openly available building density estimates, and there is insufficient coverage of these main population centres in Iran by openaccess building footprint data sets to estimate building density (Google, 2023; Microsoft, 2024). Instead, we use an estimate of population density by WorldPop (Tatem, 2017) as a proxy for residential exposure to subsidence induced hazard. We then estimate the total population exposed to each class of subsidence induced hazard. Finally, we estimate subsidence induced risk pixel‐wise for each subsiding region using the equation:

Risk = ((3 × SHG) + Vsub) × Popden (2)

where SHG is subsidence horizontal gradient or angular distortion in percent, Vsub is the absolute land subsidence velocity in m/yr, and Popden is population density (people per 100 m2 pixel) (Fernández‐Torres et al., 2020).

We additionally estimate total horizontal displacement ΔdEi and resulting horizontal strain ϵ (Cigna & Tapete, 2021b; Tandanand & Powell, 1991):

ΔdEi

ϵ = (3)

l

where ΔdEi is the change in east‐west displacement observed between two locations separated by distance l (Figure 2). Similar to our method to calculate angular distortion, we calculate horizontal strain by multiplying east‐west surface velocities by time series length and then fit a plane to each moving 3 × 3 window of pixels, again only if more than 2 pixels surround the central pixel. The value of ϵ for the central pixel is the gradient of that plane in the x‐direction divided by l. Positive ϵ indicates extension or hogging, whilst negative ϵ indicates compression or sagging (Cigna & Tapete, 2021b). We use both angular distortion and horizontal strain spatial maps and profiles to estimate the locations of previously unmapped faults.

For frames 028A_05385_191813 and 035D_05397_131013 (Tehran, Figure S4 in Supporting Information S1) we calculate angular distortion and horizontal strain using ∼50 m pixel width vertical and east‐west surface velocities, respectively, to test the impact of velocity pixel width on velocity derivative magnitudes.

2.6. Independent Component Analysis

We apply ICA to our InSAR time series to isolate elastic and inelastic temporal components of deformation within Iran's subsiding regions (Figure 2). This approach assumes that time series component sources are non‐Gaussian and statistically independent (Ebmeier, 2016; Gaddes et al., 2018), and that the temporal signatures of elastic and inelastic subsidence are consistent in time. Inelastic deformation is approximated to be linear in time as compaction is proportional to the applied, constant stress (Ezquerro et al., 2014; Haghshenas Haghighi & Motagh, 2019; Mirzadeh et al., 2023). The elastic, reversible component of deformation, is often considered seasonal, for example, in Central Valley, California (Ojha et al., 2019), but can be more complex in space and time, reaching on the order of several centimeters per year (Chaussard et al., 2014, 2017). A temporally independent component that acts in multiple subsiding regions demonstrates that a common process (such as inelastic or elastic compaction), or at least a process with the same temporal signature, may be taking place in all subsidence regions. We therefore apply temporal ICA to time series across whole LiCS frames that include several subsidence regions to test whether elastic and inelastic compaction have the same temporal signature across several, neighboring subsidence regions.

We apply a similar time series inversion methodology to generate time series for ICA to that used to produce higher resolution (100 m) InSAR time series and velocity fields. For ICA, however, we produce 1 km pixel width surface velocities and time series for entire LiCS frames (Figure S4 in Supporting Information S1; Lazecký et al., 2020) across Iran in both radar look directions. This resolution is chosen to reduce compute load and is considered sufficient to extract the dominant subsidence signals for a region. Additionally, we mask pixels where average coherence magnitude is less than 0.3 to exclude noisier velocity pixels from ICA analysis and retrieve cleaner components. We finally crop frames that span the Alborz mountain range (northern Iran, Figure 11) to only their southern extent to remove noise contaminated displacement time series and thus simplify the application of ICA. In contrast to the higher resolution subsiding region velocity estimates, for ICA we only use data acquired January 2016–December 2021. We exclude data acquired before 2016 as this period includes the Sentinel‐1A initial ramp‐up phase with inconsistent acquisition frequencies (September 2014–May 2015; e.g. Potin et al., 2015). We also exclude data acquired after the failure of Sentinel‐1B in December 2021 due to the subsequent change in Sentinel‐1 acquisition plan which reduced data acquisition frequency in Iran. After 2021, therefore, time series temporal resolution in Iran changes, which reduces the resolvability of independent components from the time series. For frames in track 159A (Figure S4 in Supporting Information S1) there is a Sentinel‐1 acquisition gap March 2017–March 2018. We limit the time series further in this track to March 2018–December 2021 to prevent this gap affecting the time series temporal pattern.

We employ temporal ICA using the FastICA algorithm (Hyvärinen, 2013; Pedregosa et al., 2011) on a frame‐byframe basis. We choose to run temporal ICA rather than spatial ICA due to the former retrieving a cleaner inelastic signal than the spatial approach, as the temporal ICA searches for consistent patterns in time rather than space. Each frame has ∼5–9 × 104 pixels at 1 km resolution with each pixel time series ∼160–200 epochs in length.

We anticipate independent processes that cause elastic and inelastic subsidence respectively to include shortterm recharge of groundwater stores by precipitation, as well as long‐term groundwater depletion. We also anticipate some independent noise sources to remain in our velocities including residual atmospheric noise and residual noise due to errors in estimating satellite acquisition location. Overall, we therefore expect at least two retrieved independent components to represent surface deformation and at least one independent noise component. We therefore conduct temporal FastICA analysis using three, four, and five components for each candidate frame.

Our approach outputs components with distinctive temporal eigenvector time series (indicating the typical temporal signature of the component) and a spatial weighting matrix (indicating how strongly weighted the temporal signature is at each pixel). After ICA, we reconstruct time series for each pixel by taking the outer product of the temporal components and corresponding mixing matrix rows (e.g., Ebmeier, 2016). Based on the reconstructed, component velocity fields and time series, we manually select the optimum number of components to use during temporal ICA for each frame. This difference in components is necessary as some frames include additional non‐subsidence deformation including earthquakes or post‐seismic relaxation (Liu et al., 2021; Watson et al., 2024). We automatically identify the inelastic component of subsidence for each frame, and manually identify the elastic component. Additional steps to identify the components are necessary as FastICA uses random initialization values to estimate the unmixing matrix (a derivative of the mixing matrix) before performing fixedpoint iteration to retrieve the independent components (Hyvärinen & Oja, 2000). This random initialization means that independent components are retrieved in a random and unpredictable order.

2.6.1. Inelastic Component Identification

To identify inelastic components automatically we assume that the inelastic component is negative with a constant rate in time. To find the component for each frame which fits this assumption, we first calculate the gradient of reconstructed time series pixel‐wise across each frame. We then measure the skew of the distribution of these time series gradients for each component. We assume and therefore identify the inelastic component as consistently having a large, negative skew to this distribution due to the dominance of a relatively small number of extremely negative pixel gradients (Table S4, Figure S7 in Supporting Information S1).

2.6.2. Elastic Component Identification

Automatic identification post ICA is more challenging for the elastic component than the inelastic component. This is firstly because the elastic signal does not always behave as a simple linear plus seasonal model. Previous work found that the elastic signal does not always have a constant annual period and can be impacted in the longer term by drought (Chaussard & Farr, 2019). Secondly, a topographically correlated atmospheric component tends to be derived from input LoS velocities when using our approach (Figure S8a–S8c in Supporting Information S1). This component has a strong seasonal trend and strongly correlates spatially with small‐scale topographic features. This signal persists despite our attempts to remove atmospheric noise using weather models via GACOS (Yu, Li, Penna, & Crippa, 2018). GACOS tends to successfully remove long wavelength atmospheric noise, with short wavelength topographically correlated noise remaining (Figure S8a–S8c in Supporting Information S1; e.g. Lohman & Simons, 2005; Ebmeier, 2016). Previous studies have demonstrated the ability of ICA to isolate atmospheric noise from InSAR time series (Belhadj‐aissa et al., 2024).

Due to these challenges, it is insufficient for automatic identification to assume that elastic deformation has an annual period. Instead, we consider components with extreme weighting matrix values in subsidence basins as candidates for elastic deformation components (e.g., Figure S9 in Supporting Information S1). Previous studies correlate components of deformation derived from ICA of InSAR time series of subsidence to groundwater well data to identify the elastic component (e.g., Chaussard & Farr, 2019; Gualandi & Liu, 2021). We could not repeat this approach in Iran due to difficulties accessing the country's groundwater well data. Data we do have access to does not provide groundwater level data at sufficient depths to analyze depletion of deep aquifers which may control elastic deformation. Moreover, the recorded groundwater levels continue only until 2019, around halfway through our ICA time series. Instead, we test the correlation of candidate elastic component deformation time series with precipitation time series (Hersbach et al., 2023), with the caveat that precipitation alone does not control groundwater levels. We divide Iran into climate zones based on Alizadeh and Keshavarz (2005) (Figures S10 and S11 in Supporting Information S1) and aggregate monthly reanalysis total precipitation data from ECMWF (Hersbach et al., 2023) for each climate zone. We assume that precipitation within each climate zone has a similar temporal pattern. We compare climate zone precipitation time series to elastic component candidate time series; where the two visually correlate we assume an element of elastic deformation is present.

For a more detailed, local analysis, we conduct watershed analysis in Tehran based on 30 m topographic data (Copernicus, 2019) to derive watersheds with a minimum area of 360,000 km2 (Figure S12 in Supporting Information S1). As we are interested in elastic deformation in the Western Tehran Plain, we select the watershed which feeds this subsidence region (Figure S5 in Supporting Information S1) for further analysis. We resample hourly reanalysis precipitation data from ECMWF (Hersbach et al., 2023) for that watershed, 2016–2021, into 12 day totals. This twelve day frequency matches the mode acquisition interval of Sentinel‐1 in Iran. We sum and normalize by area the resultant precipitation time series before examining the correlation of precipitation and candidate elastic deformation component time series.

2.6.3. Inelastic and Elastic Deformation in Space and Time

Once inelastic and elastic components are identified, we estimate how the contributions of each vary in space and time. To estimate the proportion of inelastic deformation in space, we divide the inelastic component magnitude by the LoS velocity multiplied by the time series length to calculate the inelastic/mixed signal ratio. In this scenario, any remaining, non‐inelastic signal could be attributed to residual atmospheric noise (which we assume to be small in magnitude), elastic deformation (which can range in magnitude), or deformation which is driven neither by inelastic nor elastic aquifer compaction. We compare this ratio for three main regions where we identify an elastic component: Tehran and north‐central Iran, Mashhad and north‐west Iran, and Rafsanjan and south‐central Iran (Figure 12).

To estimate the relative, proportional change of elastic to inelastic deformation over time, we first multiply elastic and inelastic eigenvector time series by the mean weighting matrix value for subsidence regions within that frame. We smooth the elastic time series by taking the mean value at the current and previous time step. We then calculate the elastic/inelastic deformation ratio (Red/id) using these reconstructed subsidence region time series (Chaussard & Farr, 2019; Ezquerro et al., 2014; Haghshenas Haghighi & Motagh, 2019).

3. Results

3.1. Subsidence Region Surface Velocities

We identify 106 subsiding regions across Iran using Sentinel‐1 InSAR time series analysis. The regions with the fastest observed vertical subsidence rates (Table S5 in Supporting Information S1, Figure 4) are located towards the centre of Iran: Rafsanjan (99th percentile vertical velocity − 221 mm/yr), Bardsir 1 (− 166 mm/yr), and Saveh (− 166 mm/yr). Twenty‐four regions experience vertical subsidence exceeding 100 mm/yr using 99th percentile values, and 34 regions if the maximum velocity is used. Qazvin is the largest subsidence region with over 2,400 km2 subsiding faster than 10 mm/yr. In total, we estimate over 31,400 km2 of Iran is subsiding (nontectonically) faster than 10 mm/yr (Table S5 in Supporting Information S1).

In general, subsidence bowls are located in valleys between steep topography (Figures 5a, 5b and 5d). Such a relationship is expected and previously observed (Herrera‐García et al., 2021) given the accumulation of sediments in valley accommodation space in which aquifers development. Indeed, the shape of subsidence regions is controlled strongly by the presence of Quaternary sediments (NIOC, 1957), with subsidence margins often matching the extent of mapped Quaternary sedimentary deposits. We additionally observe that, in some instances, previously mapped faults coincide with subsidence bowl margins (Figure 5).

Focusing on east‐west velocities (Figure 4b), only a few regions show notable horizontal deformation within areas of rapid vertical subsidence. Rafsanjan and Bardsir 1 (blue box Figures 4b; Figure 5b) is one such region. Here the west side of the subsiding region is moving eastwards, and the east side moving westwards. The change in sign of the east‐west velocities coincides with the peak in vertical velocities (Figure 5d). This pattern indicates a collapsing in of the surface towards the centre of the subsidence bowl simultaneously with vertical subsidence. The discontinuity between east and west surface velocities is not always parallel to the subsidence margin. Aside from Rafsanjan and Bardsir 1, most regions' east‐west motion is dominated by regional tectonic motion.

The dominant land cover for subsiding regions in Iran is Farm Land. We find that ∼77%, or over 24,300 km2, of land subsiding vertically faster than 10 mm/yr in Iran is covered by Farm Land (Figure 5f; Ghorbanian et al., 2020). The remaining regions are covered by kalut (a land cover type formed from soft, silty deposits through wind erosion rich in salt‐deposits, 11% Ghorbanian et al., 2020), clay (8%), urban areas (3%), and sand (0.4%). The dominance of agricultural land is consistent with several previous studies which observe rapid subsidence where groundwater is used for agricultural activities (Amighpey & Arabi, 2016; Babaee et al., 2020; Goorabi et al., 2020; Haghshenas Haghighi & Motagh, 2024; Motagh et al., 2008, 2017). Our results (Figure 5c) suggest that subsidence is faster in Farm Land covered regions (maximum rate 338 mm/yr) than those with other land cover types (Clay maximum 323 mm/yr, Urban maximum 254 mm/yr).

We observe several occurrences of salty land or kalut with apparent or real vertical uplift (e.g., Figure 6). Geologically, kalut comprises a range of deposits but includes evaporative deposits formed of lenses of gypsum, anhydrite, salty argillites, siltsones, sands, and small conglomerates (Aghanabati, 2004). Kalut therefore has a high salt content.

3.2. Comparison to Potential Subsidence Hazard

We find that ∼11% of subsidence predicted by Herrera‐García et al. (2021) is observed using Sentinel‐1 InSAR (Figure S13 in Supporting Information S1). Where there is subsidence exceeding 10 mm/yr, 92% is predicted by Herrera‐García et al. (2021). These results are conservative because we exclude apparent subsidence that does not exceed 10 mm/yr over a region larger than 4 km2. Additionally, we exclude from our subsiding regions on the southern Caspian Sea coastal areas—which Herrera‐García et al. (2021) classify as having “Very High” subsidence potential—due to our uncertainty in vertical surface velocities derived for this area (Section 2.1; Figure S2 in Supporting Information S1).

3.3. Subsidence Induced Risk

We calculate angular distortion for all 106 subsidence regions. We use this metric to calculate subsidence induced risk as well as provide evidence for inter‐ and intra‐basin geological structures that are not evidently mapped in previous work (Figure 7, Table S6 in Supporting Information S1). Across Iran, 51,500 people over 242 km2 are exposed to high levels of angular distortion.

3.3.1. Hazard and Risk in Subsiding Cities

Seven of the 10 most populated cities in Iran (Tehran, Mashhad, Isfahan, Karaj, Shiraz, Qom, Urmia) are located near one of the 106 subsiding regions.

The centre of Mashhad city is located between two subsidence bowls (dashed red lines, Figures 10a–10f), with regions of high population density (>50 people per 100 m2) extending out of the city's centre to the north‐west. In this north‐west region, vertical subsidence is at maximum 178 mm/yr. Due to steep differential subsidence gradients in this north‐west region, Mashhad has the largest area of medium (96 km2) and high (4 km2) angular distortion of any city in Iran (Table S7 in Supporting Information S1). Additionally, subsidence induced risk in north‐west Mashhad reaches ∼12 due to this area's high population density. Other cities in this study have maximum risk values of ∼25 in Tehran and Karaj; ∼8 in Isfahan; ∼5 in Shiraz; ∼2 in Qom; and ∼1.4 in Urmia.

In the country's most densely populated region, Tehran, there are two main population centres: Tehran and Karaj (Figure 1). Tehran and Karaj face some of the fastest observed vertical subsidence rates (>150 mm/yr) and highest angular distortion magnitudes (>0.085% and >0.1%, respectively) in the cities examined in this study (Figures 8(

📷

Figure 6. Example subsiding regions with adjacent uplifting areas. (a–c) Damghan; (d–f) Kavir Sirjan; (g–i) Kashan. (a, d, g) vertical velocities VU where blue indicates subsidence, red uplift. Red dashed outline highlights vertical subsidence velocities faster than 10 mm/yr. (b, e, h) land cover from Ghorbanian et al. (2020). White dashed line indicates subsidence extent. (c, f, i) 23,600 in 6.8 km2). A further 118,800 people over 59 km2 in Karaj are exposed to medium levels of this hazard. In Tehran city, around half the number of people are exposed to high surface faulting hazard as in Karaj (∼10,500 in 4.3 km2). However, due to its high population density, Tehran faces the highest observed subsidence induced risk in this study (>25), particularly to the south‐west of the city's centre (Figure 8).

Where we estimate low exposure to subsidence induced hazard and thus infer low risk to residential property, this does not exclude the exposure to non‐populated infrastructure of subsidence hazard. For example, in the densely populated centre of Qom city, angular distortion magnitudes are low (∼0.0004%), resulting in a low risk metric (∼0, Figure S15 in Supporting Information S1). However, several railway routes and roads cross the main Qom subsidence bowl to the south‐east of the city (Masoumi et al., 2022), meaning exposure and thus risk of damage to infrastructure by steep differential subsidence gradients is probably elevated. As we use population density as our exposure proxy, our risk metric does not account for non‐populated infrastructure exposure to differential subsidence hazard.

📷

Figure 7. Angular distortion mapped in space for all 106 subsiding regions. Seven of the 10 most populated cities in Iran (marked) are located near regions of subsidence. (a–f) plots angular distortion of these seven subsiding regions. Green circle in d) highlights area of noise caused by noise in the Interferometric Synthetic Aperture Radar velocities. Background Digital Elevation Model Hillshade from Esri (2024).

Shiraz and Urmia's most dense population centres are located adjacent rather than on the main subsiding region (Figures 7 and 9, Figure S16 in Supporting Information S1) posing reduced risk to the current residential population's infrastructure.

Figure 8.

3.3.2. Velocity Derivative Noise Floor

Some stable regions' angular distortion patterns appear clean and noise free, whilst other regions have a dappled pattern (green circle, Figure 8d). The distribution of noise outside of subsiding regions relates mostly to land cover and coherence. Where land cover is vegetated and coherence is low, for example in Farm Land (Figure 5), the velocity noise floor is higher meaning it is more difficult to accurately detect low values of angular distortion.

3.3.3. Impact of Pixel Width on Velocity Derivatives

We test the impact of velocity field resolution on resulting angular distortion and horizontal strain magnitudes. We find that using a finer pixel spacing (∼50 m compared to ∼100 m) for input surface velocities results in larger peak values of angular distortion, larger areas of high values of angular distortion, and thus greater estimates of people exposed to high levels of angular distortion. For example, in the Western Tehran Plain region (Figure S17c in Supporting Information S1) the ∼50 m pixel approach estimates 28,360 people in 10.7 km2 are affected by high levels of angular distortion, ∼8% more people and ∼30% greater area estimated than using ∼100 m pixel products (26,276 and 8.26 km2 respectively). This increase in angular distortion magnitude is particularly evident for pixels with already high angular distortion values such as locations of mapped faults (51.25°E, 25.55°N, Figure S17c in Supporting Information S1; faults in Figure S17d in Supporting Information S1). Using these higher resolution products does propagate a higher level of noise, such as in the mountains north of Tehran (∼0.02% difference).

3.4. Subsidence Gradient Spatial Patterns

In addition to their use in estimating subsidence induced risk, we use angular distortion (Figure 7) and horizontal strain (Figure S14 in Supporting Information S1) to identify deformation gradients that may correspond to geological structures, some of which are not evident in fault maps.

3.4.1. Angular Distortion Spatial Patterns

We use surface velocities and their derivatives in Tehran to exemplify the impact of geological structures on subsidence gradients. We also identify similar velocity and derivative patterns in other densely populated subsiding regions. In the east of the Western Tehran Plain, the pattern of high angular distortion values appears in places anastomosing, particularly to the south‐west of Tehran city (Figure 8b). Such patterns are also observed in the subsidence region north‐west of Mashhad city (Figure 10b), and north‐east of Isfahan city (Figure S18 in Supporting Information S1). In the case of Tehran's anastomosing angular distortion patterns, their strike (northnorth‐east, south‐south‐west) aligns with the strike and locations of exposed parts of the Quaternary Kahrizak formation (Mahmoudpour et al., 2016) exposed via north‐north‐east, south‐south‐west faulting within the Western Tehran Plain.

As well as anastomosing features in the Western Tehran Plain, we also observe more linear features of high (∼0.04%) angular distortion. For example, profile A‐A′, Figure 8e, spans two such linear features within the Western Tehran Plain (8b) and shows the coincidence of these peaks in angular distortion (and thus changes in vertical velocity, dark‐blue), with both a mapped fault trace (Kahrizak Fault, thick, dashed gray Figure 8e; Hessami & Jamali, 2006) and the extension of another fault trace (Eyvanekey Fault, unmapped fault splay in Figure 8e in thin dashed grey). These faults also align with a divergence in the horizontal strain pattern (profile A‐A′, Figures 8d and 8e). North and west of these mapped faults and velocity gradient excursions are several other linear features within the angular distortion pattern. These linear features have a similar strike to mapped tectonic faults and a similar value of angular distortion (∼0.06%, Figure 8b) to that where faults are mapped. In profile B‐B’ (Figure 8f), we see a similar pattern in both angular distortion and horizontal strain values to profile.

📷

Figure 8. (a) Hillshaded Copernicus 30 m Digital Elevation Model (DEM) (Copernicus, 2019) overlain with mapped faults (black; Hessami & Jamali, 2006) and outlines of subsidence faster than 10 mm/yr (red). Shaded area denotes population density exceeding 15 people per 100 m2 pixel (WorldPop, 2024). Profiles A‐A’ (d) and B‐B’

(e) are also mapped. (b) Angular distortion β for Tehran and Western Tehran Plain subsiding region. (c) Horizontal strain ϵ for the same region. (d) Sentinel‐2 5% cloud composite as in Figure 6 (e) Profiles of vertical (navy) and horizontal (red) surface velocities derived from Interferometric Synthetic Aperture Radar (this study), β (light blue) and ϵ (light red) derived from surface velocities along A‐A′ in (a–d). Grey elevation profile (A‐A′) from Copernicus 30 m DEM. Locations of a mapped fault trace (thick, grey dashed), unmapped potential extension of a mapped fault trace (thin, grey dashed), and where subsidence exceeds 10 mm/yr are also shown. (f) Same as (d) but for profile B‐B′ where no fault is mapped by Hessami and Jamali (2006). All profiled data in (e) and (f) plot the mean; first and third quartile; fifth and ninety fifth percentile of data in a 1 km wide swath along profiles A‐A′ or B‐B’.

📷

Figure 9. Vertical VU (a) and horizontal VE (e) surface velocities; angular distortion β (b); horizontal strain ϵ (f); population density (c); and risk map (d) for city Shiraz in central Iran. (a, b and c) contribute to (d) (Equation 2). Faults in (a) from Hessami and Jamali (2006). Profiled data in (g) corresponds to profile A‐A′ in (f). β (light blue line) is a derivative of VU (navy), whilst ϵ (light red) is a derivative of VE (dark red). Each line in (g) plots the mean, 5th and 95th percentile, and first and third quartile values per 1 km wide swath, spaced every 140 m, along line A‐A’.

A‐A′ at a step in topography (grey profile line). Topography here flattens where angular distortion reaches a peak of 0.06%, then resumes a similar slope to the north as in the south (dashed grey line, Figure 8f). This step corresponds to a ∼4 m vertical topographic offset (black dashed lines, Figure 8f).

We observe angular distortion linear features elsewhere in Iran. For example, profiled velocities in Figure 9 (A‐A′) reveal a sharp angular distortion peak of nearly 0.08% striking north‐west, south‐east along the western.

📷

Figure 10. Surface velocities, strain derivatives, population density, and subsequent risk map for city Mashhad in central Iran. Faults in (a) from Walker et al. (2010). Further information consistent with Figure 9.

subsidence margin (Figure 9g). This strike matches that of the surrounding topography (background, Figures 9a–9f), which is controlled by north‐west, south‐east striking anticlines to the west and east on this Quaternary aquifer (NIOC, 1957). These anticlines are composed of Eocene volcanic rocks with cores of Oligo‐Miocene and Cretaceous limestone. Similarly to other regions, the angular distortion peak corresponds with a change in sign of horizontal strain.

Within the Bardaskan subsidence bowl (Figure 15b), we observe six north‐east, south‐west trending linear angular distortion features marking zones of elevated subsidence rate. Horizontal strain in space also appears to.

📷

Figure 11. Inelastic components automatically isolated from ascending (Asc) and descending (Desc) 1 km velocity time series using Independent Component Analysis overlain on a Digital Elevation Model Hillshade from Esri (2024). (a) Inelastic components in the descending direction. Note here we plot cumulative displacement over 6 years rather than velocity. Red here is motion towards the satellite, blue is motion away. Yellow cross marks the location of pixel time series in descending frame in Figure S20 in Supporting Information S1. (b) Inelastic components in the ascending direction. Pink cross marks location of pixel time series in ascending frame in Figure S20 in Supporting Information S1. (c) Input LoS descending velocities as calculated by application of LiCSBAS (Morishita et al., 2020). Note here we plot velocities not cumulative displacement. (d) Line‐of‐sight ascending velocities calculated by LiCSBAS application.

change in magnitude at these features. To the east of this subsidence bowl near Kashmar city (Figure 15a), there are four further angular distortion features with a roughly linear pattern. On the north‐east margin of Bardaskan subsidence region there is a band of high (≥0.08%) angular distortion.

3.4.2. East‐West Velocity Spatial Patterns

Profiles of angular distortion and horizontal strain reveal in some regions (Karaj, Mashhad, Isfahan, and Qom) similar patterns of east‐west surface velocities to those observed in Rafsanjan region (Figures 4 and 5), with eastwest velocities suggesting a “collapse” of the surface towards the basin centre.

In Karaj, peak vertical subsidence rates exceed 150 mm/yr (163 mm/yr). The region has a bowl of north‐south trending rapid subsidence to the west of the city (Figure 16). Peak subsidence rates in this bowl (profile A‐A, Figure 16) coincide with a change in sign of east‐west velocities. In the subsiding region north‐west of Mashhad (Figure 10), the strike of the discontinuity between eastwards and westwards velocities is similar to the strike of the subsidence margin. Peak vertical subsidence rates in this region exceed 150 mm/yr (178 mm/yr), as in Rafsanjan and Karaj. For Isfahan, in the south‐west extension of the subsiding region (Figure S18 in Supporting Information S1), there is a sharp linear discontinuity between eastwards (red) and westwards (blue) motion that continues in a north‐east, south‐west strike towards the centre of the main subsidence region. This discontinuity is.

📷

Figure 12. Ratio of the inelastic component over the original mixed signal in the line‐of‐sight direction for ascending data in north‐central Iran (a), north‐east Iran (b), and south‐central Iran (c). These regions are the three main regions where we identify an elastic component. The ratio value indicates roughly the proportion of the signal which is due to inelastic compaction. Values exceeding one are possible due to the inelastic component having a very large, negative value that may exceed the mixed signal if the mixed signal contains other components which, during the entire time series, are positive.

📷

Figure 13. (a, d) Inelastic surface deformation component (line‐of‐sight direction, 100 m pixel width) identified automatically by our temporal Independent Component Analysis (ICA) method. Rapid deformation is limited mainly to sites of vertical subsidence faster than 10 mm/yr (pink dashed outlines). (b, e) Second “inelastic” deformation component derived by temporal ICA. Deformation is far more varied in space with both uplift and subsidence observed within subsidence bowls. Boundaries between uplift and subsidence are sometimes marked by mapped faults (black lines). (c) Time series of deformation at P1 and P2 in a and b. (f) Time series of deformation at P1 and P2 in d and e; legend the same as (c).

📷

Figure 14. (a) Line‐of‐sight displacement of a component isolated using temporal Independent Component Analysis with four components for ascending frame covering Tehran (028A_05398_191813, Figure S4 in Supporting Information S1). Red indicates motion towards the satellite, blue away. Outlines of regions of subsidence derived from 100 m velocity fields (this study) outlined in red‐dash. The restored component time series of deformation for the point (pink cross) is plotted in (b). (c) Eigenvector time series of elastic (blue) and inelastic (green) deformation for this frame at this pink point (at 1 km pixel‐width) alongside precipitation time series aggregated every 12 days to match Interferometric Synthetic Aperture Radar time series modal acquisition gap (grey, Hersbach et al., 2023).

also the location of a 0.026% peak in angular distortion. Overall, Rafsanjan, Mashhad, and Karaj peak vertical subsidence rates coincide with a change in sign of horizontal velocities, and the strike of the horizontal velocity discontinuity is often parallel to the subsidence margin or geological structures.

For Qom, the pattern of angular distortion is simpler (Figure S15 in Supporting Information S1). We observe, as identified in Rafsanjan, Isfahan, and Karaj, a discontinuity between east and west velocities (striking north‐west, south‐east) which coincides with maximum vertical velocities. A ridge of high angular distortion and the strike of surrounding topographic ridges and mapped faults (Figure S15a in Supporting Information S1; Hessami & Jamali, 2006) also follow this east‐west discontinuity.

3.5. Independent Component Analysis

We run temporal FastICA on the time series of 31 descending and 28 ascending frames across Iran (Figure 11). For 59 frames we identify an independent component with a significant negative skew to the distribution of their pixel time series gradients, which we attribute to inelastic surface deformation. This identification fails only where subsidence is masked out due to a higher coherence threshold (0.3, Figure 11). For 24 of these frames, we identify a component with some aspect of annual periodicity and extreme weighting matrix values within subsidence regions, which we attribute to elastic surface deformation within aquifers. We believe this is the first attempt to complete ICA of InSAR time series at near country scale. Figure 11 demonstrates the success of our approach in removing east‐west ramps and other noise such as earthquakes from inelastic component time series using this temporal ICA approach.

📷

Figure 15. (a) Vertical velocities VU at 100 m pixel spacing for Bardaskan region, north‐east Iran. (b) angular distortion β calculated from (a). (c) Horizontal strain ϵ calculated from east‐west velocities VE from this study at 100 m pixel width. (d) Angular distortion β derived from vertical velocities VU (a) with potential buried fault traces marked in dashed black lines. Red dashed lines in all plots mark the extent of the 10 mm/yr subsidence bowl. Straight, north‐north‐west lines of apparent high angular distortion (a, d) indicate Looking into the Continents from Space frame boundaries.

3.5.1. Inelastic Deformation

To test our success in isolating the inelastic deformation components, we compare our mixed component, observed time series to the inelastic component time series derived using temporal ICA. We select pixel time series within the Neyshabour region in north‐east Iran as the two frames which cover this region also include other, neighboring regions. In this example (Figure S20 in Supporting Information S1, left), the mixed‐signal appears to slow faster than that in the ascending direction. However, the inelastic components for both frames (right) are near linear and experience the same total cumulative displacement during the time series. As both lookdirections measure the same magnitude of inelastic displacement in the LoS, ∼100% of the direction of this inelastic deformation must be vertical.

The fastest inelastic compaction rates are located in the same subsidence region as where the fastest vertical subsidence rates are recorded. For example, Rafsanjan region in south‐central Iran records the fastest vertical subsidence rate of 342 mm/yr, 2014–2022. The same region records the greatest magnitude of inelastic compaction, 2016–2021, at ∼1,300 mm (∼216 mm/yr) in the LoS direction (Figure S21 in Supporting Information S1). This compaction magnitude is consistent between the four frames that cover this pixel, regardless of look direction. Furthermore, the subsidence region with the second largest inelastic compaction magnitude (Qazvin near Tehran), is the region with the second fastest vertical mixed signal subsidence rate. The spatial.

📷

Figure 16. Surface velocities, strain derivatives, population density, and subsequent risk maps for city Karaj in northern Iran. Faults in a) from Hessami and Jamali (2006). Further information consistent with Figure 10.

pattern of inelastic compaction within both of these subsidence bowls is also similar to the pattern of vertical subsidence. These observations are consistent with inelastic compaction driving most of the observed subsidence.

For each frame's inelastic component, we calculate the second derivative for each pixel time series of that component. We use the median pixel second derivative to compare frame‐wise changes in inelastic compaction rate, 2016–2021. We find that most mean inelastic eigenvector time series' second derivative values are small at

∼±1 × 105 eigenvector units/yr2 (Figure S22 in Supporting Information S1). The frame with the most visible change in subsidence rate also has the smallest eigenvector time series second derivative value (frame 159A_05589_131313, south‐west of Mashhad; − 0.75 x 105 eigenvector units/yr2; Figures S4 and S22 in Supporting Information S1).

We estimate the pixel‐wise proportion of total LoS signal that is inelastic compaction by calculating the inelastic/ mixed signal ratio (Figures 12a–12c) for three regions with several subsidence regions. These regions tend to have a lower inelastic to mixed signal ratio (closer to 0.6 in places) than other regions in Iran (closer to 1). Regardless of the proportion of inelastic deformation, the spatial pattern of this ratio varies within and between subsidence regions. For example, Varamin Plain (dashed box, Figure 12a) has a higher inelastic/mixed signal ratio in the west (∼0.9) than the east (∼0.5). Subsidence region margins tend to have a lower inelastic mixed signal ratio than their centres. Near Mashhad, the inelastic/mixed signal ratio varies between subsidence regions within one frame, with regions Bardaskan and Mahoolat (dashed box, Figure 12b) showing a higher proportion of inelastic deformation (∼0.9) than regions further north and east (∼0.7).

To further explore the spatial distribution of inelastic deformation components, we conduct temporal ICA on ∼100 m pixel‐width ascending LoS velocity time series for the Qazvin and Western Tehran subsiding regions (Figure S5 in Supporting Information S1). As with the 1 km product, we derive an inelastic deformation component using ICA (Figures 13a and 13d). We also identify a second deformation component for which the eignenvector time series is linear and the gradient negative. The gradient of this second linear component is shallower than the inelastic component. In space, subsidence for this second component for both Qazvin and Western Tehran Plain (Figures 13b and 13e, respectively) appears patchier, with parts of the subsiding region uplifting and parts subsiding. This component also has point‐wise, local maxima in space (Figure 13b), some of which are cone‐like with the peak of deformation relaxing radially with distance from the maximum. The boundaries between uplifting and subsiding patches within this subsidence bowl follow in some cases mapped faults (black lines). Furthermore, a sharp transition from uplift to subsidence within the Karaj subsidence region (captured in north‐west of Figure 13b) follows the strike of and is very close a ridge of high angular distortion (Figure 16).

3.5.2. Elastic Deformation

We identify an elastic deformation component for 24 frames (Figures S24 and S25; Table S4 in Supporting Information S1). Overall, most elastic component time series demonstrate near‐annual periodicity, with some changing in temporal signature around 2020.

The elastic component time series in Tehran's ascending frame closely correlates with the watershed's precipitation time series (Figure 14c, Figure S12 in Supporting Information S1). During winter months, elevated precipitation (reaching ∼5 mm/km2) is near‐synchronous with deformation towards the satellite within subsidence regions (Figure 14). Towards the end of this elastic time series (Figures 14b, 2019–2021) the proportion of elastic deformation (Red/id) decreases to 31%–36%, rather than the 41%–44% observed 2017–2018 (Figure S23 in Supporting Information S1).

4. Discussion

4.1. Nationwide Subsidence Characterisation

Our work builds on a previous Iran‐wide assessment of subsidence hazard by Haghshenas Haghighi and Motagh (2024) who used descending Sentinel‐1 interferograms and time series (2014–2020) to analyze subsidence across Iran. We significantly build on this previous work by firstly using both descending and ascending interferograms to conduct time series analysis in both look directions. We secondly decompose these time series velocities to retrieve vertical and east‐west surface velocities, 2014–2022, by incorporating north GNSS velocities. These velocity fields and subsequent velocity derivatives enable identification of subsidence gradients which we use to infer geological structures which may be previously unmapped. Finally, we use ICA to separate elastic and inelastic deformation components, and use these sources to analyze the extent of recoverable subsidence across Iran.

Compared to Haghshenas Haghighi and Motagh (2024), we estimate the area subsiding faster than 10 mm/yr in Iran is ∼21% smaller area than their reported 40,000 km2. This difference is due largely to our choice to exclude apparent observed subsidence in the Persian Gulf region of Iran. Here, similarly to regions Babol and Rasht (Figure S2 in Supporting Information S1), we observe relatively slow vertical subsidence (∼10 mm/yr) that correlates with marshland and/or shallow marine or lacustrine environments. Therefore, a large proportion of this negative velocity may be noise due to interference by vegetation or varying soil moisture. We also mask pixels with average coherence across the time series less than 0.05, which further reduces our subsidence area estimate.

We find subsidence rates have remained relatively constant between this study and previous work analyzing subsidence in Iran using SAR data. For example, Anderssohn et al. (2008) use C‐band Envisat InSAR data (2003–2006) to estimate 750–800 km2 is subsiding faster than 10 mm/yr in the Bardaskan region, with peak vertical subsidence rates 150–300 mm/yr (Figure 15). Motagh et al. (2008) use Envisat data peak subsidence reaches 270 mm/yr in this region, August–December 2004, but this estimate may incorporate a half cycle of elastic deformation. Our peak subsidence rate of 151 mm/yr is on the lower end of Anderssohn et al. (2008)'s velocity estimate. Additionally, we find the main Bardaskan subsidence bowl has grown in area by ∼40% from Anderssohn et al. (2008)'s 750–800 km2 estimate to ∼1,110 km2, whilst the smaller bowl to the south‐west adds an additional ∼160 km2 of subsidence. It is unclear whether this smaller bowl did not exist 20 years ago or it was not covered by these previous studies' data. The Yazd‐Ardakan Plain was also previously studied using Envisat data (Figure S5 in Supporting Information S1). We observe maximum vertical subsidence rates of 112 mm/yr for this region, slightly slower than the 150 mm/yr observed by Mirzadeh et al. (2021) who used Envisat and Sentinel‐1 data (2003–2020), but faster than the 94 mm/yr observed in a single Envisat interferogram by Motagh et al. (2008), 2003–2004. In the region of most rapid subsidence in Iran, Rafsanjan, we estimate peak subsidence rates of 338 mm/yr, comparable to the 300 mm/yr estimate by Motagh et al. (2017) using multiple sensors, 2004– 2016, and the 370 mm/yr by Haghshenas Haghighi and Motagh (2024).

4.2. Impact of Salt on Surface Velocities

We observe a close association of vertical uplift and mapped salty land or kalut (Figure 6). We hypothesize the mechanism of real or apparent uplift associated with salty land cover may be due to two phenomena. Firstly, salt within the kalut and salty land sediments may wet and dry cyclically over time leading to salt crystal growth which causes uplift of the ground surface without loss of phase coherence (e.g., Olivares et al., 2023). Secondly, the penetration depth of the radar signal may be influenced by soil salinity. Shao et al. (2003) experimentally show that salinity significantly affects the imaginary part of the soil's dielectric constant, especially in the 1–6 GHz frequency range, which includes the Sentinel‐1's center frequency of 5.405 GHz. They find that higher soil salinity increases the imaginary dielectric constant, leading to reduced soil penetration depth. This suggests that increased salinity might result in an apparent uplift signal in InSAR time series. However, further research is needed to confirm this hypothesis and inference, as their study primarily focuses on backscatter rather than phase changes.

4.3. Surface Velocity Derivatives to Map Geological and Structural Features

We discuss three spatial patterns of angular distortion β and horizontal strain ϵ maps (linear features, anastomosing patterns, and east‐west surface collapse) and their potential causes.

4.3.1. Linear Features

We observe linear angular distortion patterns in Western Tehran Plain, Karaj, Isfahan, Qom and Shiraz. In the Western Tehran Plain, a linear feature with ∼0.04% angular distortion follows the mapped Kahrizak Fault (Hessami & Jamali, 2006), whilst another linear feature with ∼0.06% angular distortion coincides with a potentially unmapped extension of the Eykanekey Fault (Figure 8, Profile A–A′). There are also divergences in horizontal strain at these linear features. These velocity derivative patterns suggest the mapped and potentially unmapped faults separate sediments of differing compressibilities resulting in differential subsidence rates (Amelung et al., 1999), or simply that the “unmapped” fault extension is a non‐structural boundary between two sediments of differing compressibilities. These observations support previous models by Burbey (2002, 2008), who demonstrate the association of horizontal displacements with faults in aquifer systems. Their work emphasizes the importance of monitoring horizontal displacements for information on fault locations. Our work supports their model, finding at mapped fault traces (e.g., Figure 8) a change in horizontal surface velocities and horizontal strain.

To the west of these mapped faults in Tehran, we see similar linear patterns of angular distortion and horizontal strain coincident with a ∼4 m offset in topography (Figure 8f, profile B–B′). A similar topographic offset is observed in Digital Elevation Models (DEMs) by Talebian et al. (2016) at the location the Bagh‐e‐Feyz fault trace in Tehran, which was independently located through fieldwork. Our observation of a topographic offset at the same location as a peak in angular distortion provides strong evidence for an unmapped fault. This structure may be a bedrock fault buried beneath thick Quaternary sediments that host the consolidating aquifer, possibly extending a fault trace mapped in Figure 8a. This logic could be used across Iran to identify locations of other unmapped faults within regions of rapid subsidence. Future work could include use of very high‐resolution DEMs, ICESat‐2, or other topographic data to further identify more subtle topographic fault offsets and scarps in the topography (Talebian et al., 2016) that correlate with linear features in angular distortion.

Near Shiraz, a linear peak in angular distortion follows the strike of the subsidence margin and the surrounding anticlinal topography. These observations suggest that this margin could be the boundary between compressible, Quaternary, basin infill sediments and older, incompressible, anticlinal basement rocks. Alternatively, the margin could be an aquifer‐bounding fault separating basement rocks from aquifer sediments. Field evidence may confirm the nature of this margin.

In the Bardaskan region, we observe 10 linear peaks in angular distortion, most of which strike north‐east, southwest. Such patterns align with previous observations by Anderssohn et al. (2008), who observe a complex subsidence pattern of north‐east, south‐west elongated zones of subsidence in this region (Kashmar Valley, their study). They attribute this pattern to a buried system of south‐west, north‐east striking tectonic block or graben structures beneath aquifer sediments. Anderssohn et al. (2008) mark the locations of six potential buried basement faults. We find vertical velocities (Figure 15a), angular distortion (b) and horizontal strain (c) support this buried block hypothesis and reveal angular distortion peaks in similar locations to changes in velocity magnitudes identified by Anderssohn et al. (2008). We identify an additional four peaks in angular distortion, β, outside of the original work's study region. One ≥0.08% peak on the north‐east basin margin could be caused by increasing sedimentary thickness as the basin depth increases, a change in sediment type, or a basin‐bounding fault.

4.3.2. Anastomosing Features

We observe anastomosing patterns of angular distortion that resemble channels in the Western Tehran Plain (Figure 8), Mashhad (Figure 10), and Isfahan ‐ Barkhar near Isfahan (Figure S18 in Supporting Information S1). Similar patterns have been observed in Mexico City (Cabral‐Cano et al., 2008, 2015) which coincide with transitional piedmont zones where, in this case, the geology transitions from volcanic structures to clay‐rich Quaternary lacustrine deposits. Furthermore, in Aguascalientes, Mexico, anastomosing angular distortion values higher than 0.084 correlate 50%–100% of the time with surveyed faults and fissures (Cigna & Tapete, 2021b). Other cities such as Mexico City (Cabral‐Cano et al., 2008, 2015) and Morelia, Mexico (Cigna & Tapete, 2022) also see strong correlation between high angular distortion values and mapped active faults and fissures. In the Western Tehran Plain, these anastomosing patterns align with the strike and locations of the Quaternary Kahrizak formation (Mahmoudpour et al., 2016) and related bounding faults. This formation is clayey silt, in contrast to the coarser surrounding Tehran alluvium deposit (composed of silt, sand, and pebble). This contrast in sedimentology appears therefore to induce differing subsidence rates in each formation, resulting in elevated angular distortion at their boundary. Additionally, localized formation bounding faults may affect the transmission of groundwater and contribute further to differential subsidence.

In Mashhad, similar anastomosing angular distortion patterns are observed. Due to the severe non‐linearity of these patterns, intra‐basin tectonic faults are not likely to cause these angular distortion patterns.

4.3.3. East‐West Velocities

In several locations of subsidence, east‐west velocities reveal a “collapse” of the surface towards the basin centre, such as that observed in Rafsanjan and Mashhad (Figures 4 and 5). In the subsidence bowl north‐west of Mashhad the strike of the boundary between eastward and westward velocities is similar to the strike of the subsidence margin. Additionally, the shape and peak subsidence rate (Rafsanjan 342 mm/yr, Mashhad 157 mm/yr, Table S5 in Supporting Information S1) of these subsidence bowls are similar, with both basins being elongate with high aspect ratios. These similarities suggest basin shape and subsidence rate may control the magnitude and thus observation of this east‐west surface collapse using InSAR.

We observe a similar pattern in east‐west velocities in other regions including Isfahan, Qom, Shiraz, and Karaj. In Karaj (Figure 16), peak vertical velocities coincide with a change from positive to negative horizontal velocity (Figure 16), resulting in peaks in angular distortion and horizontal strain at this location. Subsidence in Karaj is fast (maximum 163 mm/yr) supporting the hypothesis that rapid subsidence controls observation of horizontal collapse. However, unlike in Rafsanjan and Mashhad, the boundary between east and west velocities in the northern subsidence bowl in Karaj is linear and sharp, with vertical velocities showing a north‐west, south‐east trending linear feature (Figure 16). As the surface geology is consistently Quaternary sediment (NIOC, 1957; Panahi et al., 2017) across this subsidence bowl, we rule out a sedimentary control on surface velocity patterns. The linearity of the east‐west discontinuity could indicate the location of a buried, perhaps tectonic structure that controls sedimentary thickness and thus surface velocities.

Subsidence centres near Isfahan (Figure S12 in Supporting Information S1) and Qom (Figure S13 in Supporting Information S1) cities also exhibit sharp changes in horizontal surface velocity direction. In the Isfahan subsidence region, the discontinuity between eastwards and westwards motion is near‐linear. As peak subsidence rates are relatively slow in Isfahan (∼30 mm/yr), the east‐west discontinuity is linear, and there is a peak in angular distortion at this velocity discontinuity, we hypothesize that this east‐west discontinuity is structurally controlled, perhaps by a buried tectonic bedrock fault. For the Qom–Kahak region near Qom, the strike of the boundary between eastwards and westwards velocities is similar to that of the surrounding topographic ridges and mapped faults, suggesting again that there may be some buried structural control on surface velocities that relates to the geomorphology of the local region.

The ability to observe this change in direction of east‐west surface velocities could be controlled by technical limitations or artifacts in our data set, rather than being controlled by geology. For example, this horizontal motion may occur too slowly in other regions to be observed above the horizontal velocity noise threshold (up to ∼10 mm/yr in some locations). A further limitation is the lack of high‐resolution north‐south surface velocity observations for all subsidence regions. As InSAR has low sensitivity to north‐south motion, we use an interpolated, regional GNSS field to substitute and better constrain regional north‐south motion. This approach could result in residual local north‐south velocities bleeding into and interfering with local east‐west velocities, influencing the pattern of horizontal displacement.

4.4. Widespread Irreversible Subsidence

We identify inelastic components of subsidence deformation for all analyzed subsidence regions, demonstrating the widespread extent of permanent aquifer storage loss due to land subsidence across Iran. This isolation of inelastic, permanent deformation across several subsidence regions in the same InSAR frame confirms the common, irreversible mode of surface deformation and thus aquifer sediment consolidation within aquifers across Iran. We also find that inelastic deformation constitutes nearly all the observable subsidence in these locations. Where we identify an elastic component of deformation, inelastic deformation still dominates, contributing at minimum 60% of the observed deformation magnitude. For the Western Tehran Plain, Red/id never exceeds 1, again indicating the dominance of inelastic, irreversible aquifer compaction and groundwater storage loss. Several local studies of subsidence on, for example, the Yazd‐Ardakan Plain (Mirzadeh et al., 2021) and also country‐wide analysis of subsidence (Haghshenas Haghighi & Motagh, 2024) agree that inelastic compaction dominates subsidence‐associated deformation in Iran. The agreement between these works confirms the severity of groundwater storage loss in Iran. Although short‐term changes in precipitation and evapotranspiration have impacted groundwater management choices, annual mean precipitation in Iran has in fact not changed significantly since 1986 (Moshir Panahi et al., 2020; Noori et al., 2023). This apparent stability in precipitation observed at the same time as consistent inelastic subsidence widespread across Iran suggests that nationwide, unsustainable groundwater management practices are the dominant driver of declining groundwater recharge and land surface subsidence in Iran (Noori et al., 2023).

This dominance of inelastic deformation has intra‐basin variation. In the Varamin Plain subsidence region near Tehran (Figure 12), for example, the inelastic/mixed signal ratio varies from ∼0.9 in the west to ∼0.5 in the east, suggesting a higher proportion of deformation is reversible, elastic deformation in the east. The boundary between these east and west basin sections is sharp with a north‐south strike, suggesting the magnitude of this ratio has perhaps a structural and/or sedimentological control.

There is also inelastic/mixed ratio inter‐basin variation. Near Mashhad, subsidence regions Bardaskan and Mahoolat have a higher proportion of inelastic deformation (∼0.9) than regions to their north and east (∼0.7, Figure 12). This difference in space may indicate a sedimentological variation between these aquifers. For example, a higher clay to sand ratio would lead to a greater proportion of inelastic consolidation (Zhang et al., 2014).

For all regions, subsidence bowl margins tend to have a lower inelastic/mixed signal ratio than their centres. This may be due to sedimentary thickness reducing at basin margins, resulting in slower compaction (Chaussard et al., 2017). If elastic deformation magnitude remains similar at these margins to that in the centre, then the ratio will be lower.

For isolated inelastic components, at the same location in each look‐direction, displacement time series have the same displacement magnitude in the same period of time, unlike for the mixed signal (Figure S20 in Supporting Information S1). This observation supports other studies (e.g., Chaussard & Farr, 2019) in suggesting that components of deformation derived using ICA should preferably be used to derive aquifer parameters for groundwater management, rather than the observed mixed signal.

Conducting our temporal ICA approach on higher resolution (∼100 m rather than 1 km) time series data for the Western Tehran Plain reveals firstly an inelastic deformation component but also a patchy, conical pattern to the second deformation component (Figure 13b). We hypothesize this second identified deformation component relates to compaction of a different aquifer layer or layers to that or those compacting irrecoverably and represented by the components mapped in Figures 13a and 13c. Indeed, ICA has previously been used to map deformation in aquifer layers at depth (Zheng & Simons, 2023). Due to the spatial pattern of this component, we hypothesize this represents an aquifer layer which is sensitive to groundwater abstraction on a local scale, with spatial, negative point‐wise maxima reflecting local subsidence around a groundwater abstraction site. Uplifting spatial point‐wise maxima may reflect responses to ceased groundwater pumping. This same component (Figure 13b) also shows a sharp change from uplift to subsidence with the north‐east, south‐west strike of the boundary between the two regimes following the strike of the subsidence margin. High angular distortion magnitudes are also observed here (Figure 16). This ICA evidence therefore further supports the hypothesis of a buried bedrock structural control on surface velocities in this location (Section 3.4).

4.5. The Impact of Drought on Subsidence Rate

In general, we observe a correlation between elevated precipitation in winter months and uplift of the land surface (Figure 14). This correlation suggests elevated winter precipitation volumes recharge aquifer groundwater levels and lead to some uplift in the land surface which we can observe using InSAR. The lack of an apparent lag between the two time series may be due to a low clay content as a high clay composition can induce a lag between groundwater change and deformation (Zhang et al., 2014).

In the longer term, we find two pieces of evidence which suggest that subsidence rate and deformation regime may be influenced by regional climate patterns, in this case a regional drought. Firstly, as opposed to most frames which have a near constant rate of inelastic subsidence and thus very small inelastic deformation second derivative (Figure S22 in Supporting Information S1), the inelastic subsidence rate for some frames, including 159A_05589_131313, accelerates 2018–2021. Secondly, most frames' elastic component time series show annual periodicity, but some show a change in temporal signature around 2020 (including frames 028A_05390_191813, 064D_05685_131313, 093D_05476_131313, and 166D_05541_131313; Figure S24 in Supporting Information S1). In frame 028A_05398_191813 which covers Tehran (Figure 14), we quantify that the proportion of elastic deformation (Red/id) indeed decreases to 31%–36% during the drought onset (2019–2021) from a predrought value of 41%–44%, 2017–2018.

These two changes in inelastic and elastic time series may relate to the beginning in June 2020 of a severe drought that occurred across Iran and the Euphrates and Tigris basin, 2021–2023 (Otto et al., 2023). Based on the Standardized Precipitation Evapotranspiration Index, this drought is the second‐worst on record in this region. Precipitation rates fell in the short‐term during this period (Figure S11 in Supporting Information S1) which may have, according to our results, contributed to a reduction in the volume of water recharge to aquifer sediments and thus reduced the elastic controlled surface uplift in winter months (Figure 14) as well as accelerated inelastic consolidation in some regions (Figure S22 in Supporting Information S1). This decrease in Red/id during the drought indicates an increasing dominance of inelastic deformation and reduced ability of aquifer system storage to recover elastically. However, as both during and pre‐drought Red/id is less than one, inelastic deformation is still persistently dominant throughout the observation period. The reduced elastic deformation proportion near Tehran during the drought aligns with the observations of Chaussard and Farr (2019) who, during a drought period in San Joaquin Valley, 2015–2019, observed a reduced proportion of elastic deformation. Future work should focus on characterisation of the properties of these aquifers in Iran to confirm if changes in precipitation rates are able to influence both aquifer recharge and surface deformation magnitudes on a seasonal timescale (Murray & Lohman, 2018; Shirzaei et al., 2019). Additionally, the relative impact of snow and rainfall recharge (both included in precipitation in this study) on Red/id should be investigated. Previous work by Safarianzengir et al. (2020), for example, noted a high correlation between decreasing snow fall rates and lowering groundwater levels in Iran since 2008.

The inelastic component appears to show a minor amount of seasonal signal (Figure 14c), which may reflect either imperfectly retrieved elastic signal remaining in this component (“ghosting”), or a small (1%–5%), real, recoverable component of inelastic deformation that can be realized when stresses are reduced (Sneed & Galloway, 2000).

We note that our temporal ICA approach does not identify an elastic component for several frames, precluding their ability to be assessed for drought impact. This may be because either the magnitude of the near‐annual period signal is smaller or similar to the noise threshold, or there is no real elastic signal to the subsidence deformation. The latter option might be due to a lower ratio of aquitards to aquifers leading to a smaller magnitude elastic component (Chaussard & Farr, 2019).

5. Conclusion

Sentinel‐1 InSAR time series data for 106 subsiding regions across Iran in this study provides for the first time an 8‐year (2014–2022) history of surface deformation in two component directions (ascending and descending). For each region, we derive vertical and east‐west surface velocity fields at 100 m pixel width from these InSAR time series data. We find 24 regions experience vertical subsidence rates exceeding 100 mm/yr, using 99th percentile values. As a whole, this data set provides a valuable tool for assessing the current locations of subsidence and how subsidence rates have changed over time. Through calculation of angular distortion and horizontal strain using vertical and east‐west component velocities, respectively, the locations of elevated differential settlement hazard are mapped for each subsiding region. Across the country, more than 51,500 people over 242 km2 are exposed to high levels of angular distortion. We demonstrate the utility of these surface velocity spatial derivatives in potentially identifying geological structures which are unapparent in existing maps. ICA of InSAR velocity time series is conducted on a country‐wide scale for the first time to automatically and manually isolate inelastic and elastic components of subsidence deformation, respectively. We demonstrate the dominance of inelastic deformation over elastic, highlighting the widespread aquifer storage loss across the country. We additionally demonstrate the response of elastic and inelastic deformation to changing precipitation patterns during a regional drought. We recommend future work focuses on further constraint of hypothesized, unmapped fault locations identified in this study using, for example, satellite topographic data and fieldwork investigations. Finally, we recommend characterization of aquifer properties to determine if aquifer recharge and subsequent surface deformation is impacted by precipitation changes on a seasonal timescale. The systematic application of these techniques on a regional‐scale highlights the value and potential of regional InSAR for informing nationwide and local groundwater management practices.

Reply to this discussion

Conservation

Transmission

Recovery

Visual

Accelerators

Following

Share

Top of Form

References:

Aghanabati, A. (2004). The geology of Iran. Iran’s Geology Institution Press.

Alizadeh, A., & Keshavarz, A. (2005). Status of agricultural water use in Iran, in water conservation, reuse, and recycling: Proceedings of an Iranian‐American workshop (Vol. 4, pp. 94–105). National Academies Press.

Amelung, F., Galloway, D. L., Bell, J. W., Zebker, H. A., & Laczniak, R. J. (1999). Sensing the ups and downs of Las Vegas: InSAR reveals structural control of land subsidence and aquifer‐system deformation. Geology, 27(6), 483–486. https://doi.org/10.1130/0091‐7613(1999) 027%C2%A10483:STUADO%C2%BF2.3.CO;2

Amighpey, M., & Arabi, S. (2016). Studying land subsidence in Yazd province, Iran, by integration of InSAR and levelling measurements. Remote Sensing Applications: Society and Environment, 4, 1–8. https://doi.org/10.1016/j.rsase.2016.04.001

Anderssohn, J., Wetzel, H.‐U., Walter, T. R., Motagh, M., Djamour, Y., & Kaufmann, H. (2008). Land subsidence pattern controlled by old alpine basement faults in the Kashmar valley, northeast Iran: Results from InSAR and levelling. Geophysical Journal International, 174(1), 287–294. https://doi.org/10.1111/j.1365‐246X.2008.03805.x

Ashraf, S., Nazemi, A., & AghaKouchak, A. (2021). Anthropogenic drought dominates groundwater depletion in Iran. Scientific Reports, 11(1), 9135. https://doi.org/10.1038/s41598‐021‐88522‐y

Babaee, S., Mousavi, Z., Masoumi, Z., Malekshah, A. H., Roostaei, M., & Aflaki, M. (2020). Land subsidence from interferometric SAR and groundwater patterns in the Qazvin plain, Iran. International Journal of Remote Sensing, 41(12), 4780–4798. https://doi.org/10.1080/ 01431161.2020.1724345

Baldi, P., Casula, G., Cenni, N., Loddo, F., & Pesci, A. (2009). GPS‐Based monitoring of land subsidence in the Po plain (northern Italy). Earth and Planetary Science Letters, 288(1), 204–212. https://doi.org/10.1016/j.epsl.2009.09.023

Belhadj‐aissa, S., Simard, M., Jones, C. E., Oliver‐Cabrera, T., & Christensen, A. (2024). Separation of water level change from atmospheric artifacts through application of independent component analysis to InSAR time series. Earth and Space Science, 11(7), e2024EA003540. https://doi.org/10.1029/2024EA003540

Biggs, J., Wright, T., Lu, Z., & Parsons, B. (2007). Multi‐interferogram method for measuring interseismic deformation: Denali fault, Alaska. Geophysical Journal International, 170(3), 1165–1179. https://doi.org/10.1111/j.1365‐246X.2007.03415.x

Bouwer, H. (1977). Land subsidence and cracking due to ground‐water depletiona. Groundwater Series, 15(5), 358–364. https://doi.org/10.1111/j. 1745‐6584.1977.tb03180.x

Burbey, T. J. (2002). The influence of faults in basin‐fill deposits on land subsidence, Las Vegas valley, Nevada, USA. Hydrogeology Journal, 10(5), 525–538. https://doi.org/10.1007/s10040‐002‐0215‐7

Burbey, T. J. (2008). The influence of geologic structures on deformation due to ground water withdrawal. Groundwater Series, 46(2), 202–211. https://doi.org/10.1111/j.1745‐6584.2007.00395.x

Cabral‐Cano, E., Dixon, T. H., Miralles‐Wilhelm, F., Diaz‐Molina, O., Sanchez‐Zamora, O., & Carande, R. E. (2008). Space geodetic imaging of rapid ground subsidence in Mexico City. Geological Society of America Bulletin, 120(11–12), 1556–1566. https://doi.org/10.1130/B26001.1

Cabral‐Cano, E., Solano‐Rojas, D., Oliver‐Cabrera, T., Wdowinski, S., Chaussard, E., Salazar‐Tlaczani, L., et al. (2015). Satellite geodesy tools for ground subsidence and associated shallow faulting hazard assessment in central Mexico, In Proceedings of the International Association of Hydrological Sciences, Proceedings of IAHS, Vol. 372, pp. 255–260. Copernicus GmbH. https://doi.org/10.5194/piahs‐372‐255‐2015

Castellazzi, P., Garfias, J., Martel, R., Brouard, C., & Rivera, A. (2017). InSAR to support sustainable urbanization over compacting aquifers: The case of Toluca valley, Mexico. International Journal of Applied Earth Observation and Geoinformation, 63, 33–44. https://doi.org/10.1016/j. jag.2017.06.011

Chaussard, E., Amelung, F., Abidin, H., & Hong, S.‐H. (2013). Sinking cities in Indonesia: ALOS PALSAR detects rapid subsidence due to groundwater and gas extraction. Remote Sensing of Environment, 128, 150–161. https://doi.org/10.1016/j.rse.2012.10.015

Chaussard, E., & Farr, T. G. (2019). A new method for isolating elastic from inelastic deformation in aquifer systems: Application to the San Joaquin valley, CA. Geophysical Research Letters, 46(19), 10800–10809. https://doi.org/10.1029/2019GL084418

Chaussard, E., Milillo, P., Bürgmann, R., Perissin, D., Fielding, E. J., & Baker, B. (2017). Remote sensing of ground deformation for monitoring groundwater management practices: Application to the Santa Clara valley during the 2012–2015 California drought. Journal of Geophysical Research: Solid Earth, 122(10), 8566–8582. https://doi.org/10.1002/2017JB014676

Chaussard, E., Wdowinski, S., Cabral‐Cano, E., & Amelung, F. (2014). Land subsidence in central Mexico detected by ALOS InSAR time‐series. Remote Sensing of Environment, 140, 94–106. https://doi.org/10.1016/j.rse.2013.08.038

Cigna, F., & Tapete, D. (2021). Present‐day land subsidence rates, surface faulting hazard and risk in Mexico City with 2014–2020 Sentinel‐1 IW InSAR. Remote Sensing of Environment, 253, 112–161. https://doi.org/10.1016/j.rse.2020.112161

Cigna, F., & Tapete, D. (2021). Satellite InSAR survey of structurally‐controlled land subsidence due to groundwater exploitation in the Aguascalientes valley, Mexico. Remote Sensing of Environment, 254, 112–254. https://doi.org/10.1016/j.rse.2020.112254

Cigna, F., & Tapete, D. (2022). Urban growth and land subsidence: Multi‐decadal investigation using human settlement data and satellite InSAR in Morelia, Mexico. Science of The Total Environment, 811, 152–211. https://doi.org/10.1016/j.scitotenv.2021.152211

Copernicus. (2019). Copernicus DEM ‐ Global and European digital elevation model | Copernicus data space ecosystem. https://doi.org/10.5270/ ESA‐c5d3d65

De Zan, F., Zonno, M., & López‐Dekker, P. (2015). Phase inconsistencies and multiple scattering in SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing, 53(12), 6608–6616. https://doi.org/10.1109/TGRS.2015.2444431

Ebmeier, S. K. (2016). Application of independent component analysis to multitemporal InSAR data with volcanic case studies. Journal of Geophysical Research: Solid Earth, 121(12), 8970–8986. https://doi.org/10.1002/2016JB013765

English, P. W. (1968). The origin and spread of qanats in the old world. Proceedings of the American Philosophical Society. (Vol. 112(3), pp. 170– 181). American Philosophical Society

ESA. (2021). Copernicus Sentinel‐2 (processed by ESA). https://doi.org/10.5270/S2_‐742ikth

Esri. (2024). World hillshade. [Map] Retrieved from https://services.arcgisonline.com/arcgis/rest/services/Elevation/World_hillshade/ MapServer

Ezquerro, P., Herrera, G., Marchamalo, M., Tomás, R., Béjar‐Pizarro, M., & Martínez, R. (2014). A quasi‐elastic aquifer deformational behavior: Madrid aquifer case study. Journal of Hydrology, 519, 1192–1204. https://doi.org/10.1016/j.jhydrol.2014.08.040

Famiglietti, J. S., & Ferguson, G. (2021). The hidden crisis beneath our feet. Science, 372(6540), 344–345. https://doi.org/10.1126/science. abh2867

Fernández‐Torres, E., Cabral‐Cano, E., Solano‐Rojas, D., Havazli, E., & Salazar‐Tlaczani, L. (2020). Land subsidence risk maps and InSAR based angular distortion structural vulnerability assessment: An example in Mexico City, In Proceedings of the International Association of Hydrological Sciences, Proceedings of IAHS, Vol. 382, pp. 583–587. Copernicus GmbH. https://doi.org/10.5194/piahs‐382‐583‐2020

Fernández‐Torres, E. A., Cabral‐Cano, E., Novelo‐Casanova, D. A., Solano‐Rojas, D., Havazli, E., & Salazar‐Tlaczani, L. (2022). Risk assessment of land subsidence and associated faulting in Mexico City using InSAR. Natural Hazards, 112(1), 37–55. https://doi.org/10.1007/ s11069‐021‐05171‐0

Foroughnia, F., Nemati, S., Maghsoudi, Y., & Perissin, D. (2019). An iterative PS‐InSAR method for the analysis of large spatio‐temporal baseline data stacks for land subsidence estimation. International Journal of Applied Earth Observation and Geoinformation, 74, 248–258. https://doi.org/10.1016/j.jag.2018.09.018

Gaddes, M. E., Hooper, A., Bagnardi, M., Inman, H., & Albino, F. (2018). Blind signal separation methods for InSAR: The potential to automatically detect and monitor signals of volcanic deformation. Journal of Geophysical Research: Solid Earth, 123(11), 10226–10251. https:// doi.org/10.1029/2018JB016210

Galloway, D. L., Hudnut, K. W., Ingebritsen, S. E., Phillips, S. P., Peltzer, G., Rogez, F., & Rosen, P. A. (1998). Detection of aquifer system compaction and land subsidence using interferometric synthetic aperture radar, antelope valley, Mojave desert, California. Water Resources Research, 34(10), 2573–2585. https://doi.org/10.1029/98WR01285

Ghorbanian, A., Kakooei, M., Amani, M., Mahdavi, S., Mohammadzadeh, A., & Hasanlou, M. (2020). Improved land cover map of Iran using sentinel imagery within google Earth engine and a novel automatic workflow for land cover classification using migrated training samples.

ISPRS Journal of Photogrammetry and Remote Sensing, 167, 276–288. https://doi.org/10.1016/j.isprsjprs.2020.07.013

Google. (2023). Google open buildings, version 3. Retrieved from https://sites.research.google/gr/open‐buildings/

Goorabi, A., Karimi, M., Yamani, M., & Perissin, D. (2020). Land subsidence in Isfahan metropolitan and its relationship with geological and geomorphological settings revealed by Sentinel‐1A InSAR observations. Journal of Arid Environments, 181, 104–238. https://doi.org/10.1016/ j.jaridenv.2020.104238

Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D., & Moore, R. (2017). Google Earth engine: Planetary‐scale geospatial analysis for everyone. Remote Sensing of Environment, 202, 18–27. https://doi.org/10.1016/j.rse.2017.06.031

Gualandi, A., & Liu, Z. (2021). Variational Bayesian independent component analysis for InSAR displacement time‐series with application to central California, USA. Journal of Geophysical Research: Solid Earth, 126(4), e2020JB020845. https://doi.org/10.1029/2020JB020845

Haghshenas Haghighi, M., & Motagh, M. (2019). Ground surface response to continuous compaction of aquifer system in Tehran, Iran: Results from a long‐term multi‐sensor InSAR analysis. Remote Sensing of Environment, 221, 534–550. https://doi.org/10.1016/j.rse.2018.11.003

Haghshenas Haghighi, M., & Motagh, M. (2024). Uncovering the impacts of depleting aquifers: A remote sensing analysis of land subsidence in Iran. Science Advances, 10(19), eadk3039. https://doi.org/10.1126/sciadv.adk3039

Herrera‐García, G., Ezquerro, P., Tomás, R., Béjar‐Pizarro, M., López‐Vinielles, J., Rossi, M., et al. (2021). Mapping the global threat of land subsidence. Science (New York, N.Y.), 371(6524), 34–36. https://doi.org/10.1126/science.abb8549

Herring, T. (1999). Geodetic applications of GPS. Proceedings of the IEEE, 87(1), 92–110. https://doi.org/10.1109/5.736344

Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horanyi, A., Muñoz Sabater, J., et al. (2023). ERA5 monthly averaged data on single levels from 1940 to present. https://doi.org/10.24381/CDS.F17050D7

Hessami, K., & Jamali, F. (2006). Explanatory notes to the map of major active faults of Iran. Journal of Seismology and Earthquake Engineering, 8(1), 1–11.

Hooper, A., Wright, T. J., Spaans, K., Elliott, J., Weiss, J. R., Bagnardi, M., et al. (2018). Global monitoring of fault zones and volcanoes with Sentinel‐1 [Dataset]. IGARSS 2018 ‐ 2018 IEEE International Geoscience and Remote Sensing Symposium, 1566–1568. https://doi.org/10. 1109/IGARSS.2018.8517519

Hyvärinen, A. (2013). Independent component analysis: Recent advances. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(1984), 20110534. https://doi.org/10.1098/rsta.2011.0534

Hyvärinen, A., & Oja, E. (2000). Independent component analysis: Algorithms and applications. Neural Networks, 13(4), 411–430. https://doi. org/10.1016/S0893‐6080(00)00026‐5

Jasechko, S., Seybold, H., Perrone, D., Fan, Y., Shamsudduha, M., Taylor, R. G., et al. (2024). Rapid groundwater decline and some cases of recovery in aquifers globally. Nature, 625(7996), 715–721. https://doi.org/10.1038/s41586‐023‐06879‐8

Khorrami, F., Vernant, P., Masson, F., Nilfouroushan, F., Mousavi, Z., Nankali, H., et al. (2019). An up‐to‐date crustal deformation map of Iran using integrated campaign‐mode and permanent GPS velocities. Geophysical Journal International, 217(2), 832–843. https://doi.org/10.1093/ gji/ggz045

Lazecký, M., Ou, Q., Shen, L., McGrath, J., Payne, J., Espín, P., et al. (2024). Strategies for improving and correcting unwrapped interferograms implemented in LiCSBAS. Procedia Computer Science, 239, 2408–2412. https://doi.org/10.1016/j.procs.2024.06.435

Lazecký, M., Spaans, K., González, P. J., Maghsoudi, Y., Morishita, Y., Albino, F., et al. (2020). LiCSAR: An automatic InSAR tool for measuring and monitoring tectonic and volcanic activity. Remote Sensing, 12(15), 2430. https://doi.org/10.3390/rs12152430

Leake, S. A. (1990). Interbed storage changes and compaction in models of regional groundwater flow. Water Resources Research, 26(9), 1939– 1950. https://doi.org/10.1029/WR026i009p01939

Liu, F., Elliott, J. R., Craig, T. J., Hooper, A., & Wright, T. J. (2021). Improving the resolving power of InSAR for earthquakes using time series: A case study in Iran. Geophysical Research Letters, 48(14), e2021GL093043. https://doi.org/10.1029/2021GL093043

Liu, Y., Li, J., Fang, Z. N., Rashvand, M., & Griffin, T. (2020). The secondary consolidation (creep) due to geohistorical overburden pressure in the Houston‐Galveston region, Texas. In Proceedings of the International Association of Hydrological Sciences, Proceedings of IAHS, Vol. 382, pp. 315–320. Copernicus GmbH. https://doi.org/10.5194/piahs‐382‐315‐2020

Lohman, R. B., & Simons, M. (2005). Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling. Geochemistry, Geophysics, Geosystems, 6(1), Q01007. https://doi.org/10.1029/2004GC000841

Madani, K., AghaKouchak, A., & Mirchi, A. (2016). Iran’s socio‐economic drought: Challenges of a water‐bankrupt nation. Iranian Studies, 49(6), 997–1016. https://doi.org/10.1080/00210862.2016.1259286

Maghrebi, M., Noori, R., Partani, S., Araghi, A., Barati, R., Farnoush, H., & Torabi Haghighi, A. (2021). Iran’s groundwater hydrochemistry. Earth and Space Science, 8(8), e2021EA001793. https://doi.org/10.1029/2021EA001793

Maghsoudi, Y., Hooper, A. J., Wright, T. J., Lazecky, M., & Ansari, H. (2022). Characterizing and correcting phase biases in short‐term, multilooked interferograms. Remote Sensing of Environment, 275, 113022. https://doi.org/10.1016/j.rse.2022.113022

Mahmoudpour, M., Khamehchiyan, M., Nikudel, M., & Gassemi, M. (2013). Characterization of regional land subsidence induced by groundwater withdrawals in Tehran, Iran. Geopersia, 3(2), 49–62. https://doi.org/10.22059/jgeope.2013.36014

Mahmoudpour, M., Khamehchiyan, M., Nikudel, M. R., & Ghassemi, M. R. (2016). Numerical simulation and prediction of regional land subsidence caused by groundwater exploitation in the southwest plain of Tehran, Iran. Engineering Geology, 201, 6–28. https://doi.org/10. 1016/j.enggeo.2015.12.004

Mashayekhi, A. (2016). Tehran, the scene of modernity in the pahlavi dynasty: Modernisation and urbanisation processes 1925–1979. In F. F. Arefian & S. H. I. Moeini (Eds.), Urban change in Iran: Stories of rooted histories and ever‐accelerating developments (pp. 103–119). Springer International Publishing.

Masoumi, Z., Mousavi, Z., & Hajeb, Z. (2022). Long‐term investigation of subsidence rate and its environmental effects using the InSAR technique and geospatial analyses. Geocarto International, 37(24), 7161–7185. https://doi.org/10.1080/10106049.2021.1964616

Microsoft. (2024). GlobalMLBuildingFootprints. Retrieved from https://planetarycomputer.microsoft.com/dataset/ms‐buildings

Mirzadeh, S. M. J., Jin, S., Chaussard, E., Bürgmann, R., Rezaei, A., Ghotbi, S., & Braun, A. (2023). Transition and drivers of elastic to inelastic deformation in the Abarkuh plain from InSAR multi‐sensor time series and hydrogeological data. Journal of Geophysical Research: Solid Earth, 128(7), e2023JB026430. https://doi.org/10.1029/2023JB026430

Mirzadeh, S. M. J., Jin, S., Parizi, E., Chaussard, E., Bürgmann, R., Delgado Blasco, J. M., et al. (2021). Characterization of irreversible land subsidence in the Yazd‐Ardakan plain, Iran from 2003 to 2020 InSAR time series. Journal of Geophysical Research: Solid Earth, 126(11), e2021JB022258. https://doi.org/10.1029/2021JB022258

Morishita, Y., Lazecky, M., Wright, T. J., Weiss, J. R., Elliott, J. R., & Hooper, A. (2020). LiCSBAS: An open‐source InSAR time series analysis package integrated with the LiCSAR automated Sentinel‐1 InSAR processor. Remote Sensing, 12(3), 424. https://doi.org/10.3390/rs12030424

Moshir Panahi, D., Kalantari, Z., Ghajarnia, N., Seifollahi‐Aghmiuni, S., & Destouni, G. (2020). Variability and change in the hydro‐climate and water resources of Iran over a recent 30‐year period. Scientific Reports, 10(1), 7450. https://doi.org/10.1038/s41598‐020‐64089‐y

Motagh, M., Shamshiri, R., Haghshenas Haghighi, M., Wetzel, H.‐U., Akbari, B., Nahavandchi, H., et al. (2017). Quantifying groundwater exploitation induced subsidence in the Rafsanjani plain, southeastern Iran, using InSAR time‐series and in situ measurements. Engineering Geology, 218, 134–151. https://doi.org/10.1016/j.enggeo.2017.01.011

Motagh, M., Walter, T. R., Sharifi, M. A., Fielding, E., Schenk, A., Anderssohn, J., & Zschau, J. (2008). Land subsidence in Iran caused by widespread water reservoir overexploitation. Geophysical Research Letters, 35(16), L16403. https://doi.org/10.1029/2008GL033814

Müller Schmied, H., Cáceres, D., Eisner, S., Flörke, M., Herbert, C., Niemann, C., et al. (2021). The global water resources and use model WaterGAP v2.2d: Model description and evaluation. Geoscientific Model Development, 14(2), 1037–1079. https://doi.org/10.5194/gmd‐141037‐2021

Murray, K. D., & Lohman, R. B. (2018). Short‐lived pause in central California subsidence after heavy winter precipitation of 2017. Science Advances, 4(8), eaar8144. https://doi.org/10.1126/sciadv.aar8144

NIOC. (1957). Geological map of Iran. ‐ ESDAC ‐ European commission. Retrieved from https://esdac.jrc.ec.europa.eu/content/geological‐mapiran

Noori, R., Maghrebi, M., Jessen, S., Bateni, S. M., Heggy, E., Javadi, S., et al. (2023). Decline in Iran’s groundwater recharge. Nature Communications, 14(1), 6674. https://doi.org/10.1038/s41467‐023‐42411‐2

Noori, R., Maghrebi, M., Mirchi, A., Tang, Q., Bhattarai, R., Sadegh, M., et al. (2021). Anthropogenic depletion of Iran’s aquifers. Proceedings of the National Academy of Sciences, 118(25), e2024221118. https://doi.org/10.1073/pnas.2024221118

Ojha, C., Werth, S., & Shirzaei, M. (2019). Groundwater loss and aquifer system compaction in san joaquin valley during 2012–2015 drought. Journal of Geophysical Research: Solid Earth, 124(3), 3127–3143. https://doi.org/10.1029/2018JB016083

Olivares, L., Jordan, T. E., Philpot, W. D., & Lohman, R. B. (2023). Recognition of whole‐landscape changes due to extreme rain events in a hyperarid desert. Remote Sensing Applications: Society and Environment, 29, 100–927. https://doi.org/10.1016/j.rsase.2023.100927

Otto, F., Clarke, B., Rahimi, M., Zachariah, M., Barnes, C., Kimutai, J., et al. (2023). Human‐induced climate change compounded by socioeconomic water stressors increased severity of drought in Syria, Iraq and Iran. World Weather Attribution Scientific Report. https://doi.org/ 10.25561/107370

Ou, Q., Daout, S., Weiss, J. R., Shen, L., Lazecký, M., Wright, T. J., & Parsons, B. E. (2022). Large‐scale interseismic strain mapping of the NE Tibetan Plateau from Sentinel‐1 interferometry. Journal of Geophysical Research: Solid Earth, 127(6), e2022JB024176. https://doi.org/10. 1029/2022JB024176

Panahi, M. R., Mousavi, S. M., & Rahimzadegan, M. (2017). Delineation of groundwater potential zones using remote sensing, GIS, and AHP technique in Tehran–Karaj plain, Iran. Environmental Earth Sciences, 76(23), 792. https://doi.org/10.1007/s12665‐017‐7126‐3

Payne, J., Watson, A., Maghsoudi, Y., Ebmeier, S., Rigby, R., Lazecky, M., et al. (2024). Accompanying software for pre‐print: ‘widespread extent of irrecoverable aquifer depletion revealed by country‐wide analysis of land surface subsidence hazard in Iran [Software]. https://doi.org/ 10.5281/zenodo.14900087

Payne, J., Watson, A., Maghsoudi, Y., Ebmeier, S. K., Rigby, R., Lazecký, M., et al. (2024). Accompanying data for pre‐print: “widespread extent of irrecoverable aquifer depletion revealed by country‐wide analysis of land surface subsidence hazard in Iran” [Dataset]. Zenodo. https://doi. org/10.5281/zenodo.13754200

Payne, J., Watson, A., Maghsoudi, Y., Lazecký, M., Watson, S., & Elliott, J. R. (2023). Nationwide assessment of subsidence induced hazard in Iran using Sentinel‐1 InSAR [Dataset]. AGU Fall Meeting Abstracts, 2023. https://comet‐subsidencedb.org/

Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., et al. (2011). Scikit‐learn: Machine learning in python. Journal of Machine Learning Research, 12(85), 2825–2830.

Potin, P., Rosich, P., Miranda, N., Grimont, P., Bargellini, P., Monjoux, E., et al. (2015). Sentinel‐1 mission status. In 2015 IEEE international geoscience and remote sensing symposium (IGARSS) (pp. 2820–2823). https://doi.org/10.1109/IGARSS.2015.7326401

QGIS.org (2024). QGIS geographic information system. QGIS Association. Retrieved from https://www.qgis.org

Safarianzengir, V., Mahmoudi, L., Meresht, R. M., Abad, B., Rajabi, K., & Kianian, M. (2020). Monitoring and analysis of changes in the depth and surface area snow of the Mountains in Iran using remote sensing data. Journal of the Indian Society of Remote Sensing, 48(11), 1479–1494. https://doi.org/10.1007/s12524‐020‐01145‐0

Sanabria, M. P., Guardiola‐Albert, C., Tomás, R., Herrera, G., Prieto, A., Sánchez, H., & Tessitore, S. (2014). Subsidence activity maps derived from DInSAR data: Orihuela case study. Natural Hazards and Earth System Sciences, 14(5), 1341–1360. https://doi.org/10.5194/nhess‐141341‐2014

SCI. (2016). Data and statistical information. Retrieved from https://www.amar.org.ir/statistical‐information

Shao, Y., Hu, Q., Guo, H., Lu, Y., Dong, Q., & Han, C. (2003). Effect of dielectric properties of moist salinized soils on backscattering coefficients extracted from RADARSAT image. IEEE Transactions on Geoscience and Remote Sensing, 41(8), 1879–1888. https://doi.org/10.1109/TGRS. 2003.813499

Shirzaei, M., Ojha, C., Werth, S., Carlson, G., & Vivoni, E. R. (2019). Comment on “Short‐lived pause in Central California subsidence after heavy winter precipitation of 2017”. In K. D. Murray & R. B. Lohman (Eds.). Science Advances, Vol. 5(6), eaav8038. https://doi.org/10.1126/ sciadv.aav8038

Skempton, A. W., & Macdonald, D. H. (1956). The allowable settlements of buildings. Proceedings ‐ Institution of Civil Engineers, 5(6), 727– 768. https://doi.org/10.1680/ipeds.1956.12202

Sneed, M., & Galloway, D. (2000). Aquifer‐system compaction and land subsidence: Measurements, analyses, and Simulations—The holly site, Edwards air force base. U.S. Geological Survey.

Sokhadze, G., Floyd, M., Godoladze, T., King, R., Cowgill, E. S., Javakhishvili, Z., et al. (2018). Active convergence between the lesser and greater Caucasus in Georgia: Constraints on the tectonic evolution of the lesser–greater Caucasus Continental collision. Earth and Planetary Science Letters, 481, 154–161. https://doi.org/10.1016/j.epsl.2017.10.007

Solano Rojas, D. E., Wdowinski, S., Cabral‐Cano, E., & Osmanoglu, B. (2017). Differential subsidence in Mexico City and implications to its collective transport system (metro). AGU Fall Meeting Abstracts, 2017, NH32A–08.

Stephenson, O. L., Liu, Y.‐K., Yunjun, Z., Simons, M., Rosen, P., & Xu, X. (2022). The impact of plate motions on long‐wavelength InSARDerived velocity fields. Geophysical Research Letters, 49(21), e2022GL099835. https://doi.org/10.1029/2022GL099835

Susilo, S., Salman, R., Hermawan, W., Widyaningrum, R., Wibowo, S. T., Lumban‐Gaol, Y. A., et al. (2023). GNSS land subsidence observations along the northern coastline of java, Indonesia. Scientific Data, 10(1), 421. https://doi.org/10.1038/s41597‐023‐02274‐0

Talebian, M., Copley, A., Fattahi, M., Ghorashi, M., Jackson, J., Nazari, H., et al. (2016). Active faulting within a megacity: The geometry and slip rate of the pardisan thrust in central Tehran, Iran. Geophysical Journal International, 207(3), 1688–1699. https://doi.org/10.1093/gji/ggw347

Tandanand, S., & Powell, L. R. (1991). Determining horizontal displacement and strains due to subsidence. Tech. rep., Bureau of the Mines, United States Department of the Interior, United States.

Tatem, A. J. (2017). WorldPop, open data for spatial demography. Scientific Data, 4(1), 170004. https://doi.org/10.1038/sdata.2017.4 Terzaghi, K. (1925). Principles of soil mechanics, IV ‐ Settlement and consolidation of clay. Engineering News‐Record, 95(3), 874–878.

Walker, R. T., Talebian, M., Saiffori, S., Sloan, R. A., Rasheedi, A., MacBean, N., & Ghassemi, A. (2010). Active faulting, earthquakes, and restraining Bend development near Kerman city in southeastern Iran. Journal of Structural Geology, 32(8), 1046–1060. https://doi.org/10. 1016/j.jsg.2010.06.012

Watson, A. R., Elliott, J. R., Lazecký, M., Maghsoudi, Y., McGrath, J. D., & Walters, R. J. (2024). An InSAR‐GNSS velocity field for Iran. Geophysical Research Letters, 51(10), e2024GL108440. https://doi.org/10.1029/2024GL108440

Watson, A. R., Elliott, J. R., & Walters, R. J. (2022). Interseismic strain accumulation across the main recent fault, SW Iran, from Sentinel‐1 InSAR observations. Journal of Geophysical Research: Solid Earth, 127(2), e2021JB022674. https://doi.org/10.1029/2021JB022674

WorldPop. (2024). WorldPop global project population data: Estimated residential population per 100x100m grid square | Earth engine data catalog | google for developers.

Yu, C., Li, Z., & Penna, N. T. (2018). Interferometric synthetic aperture radar atmospheric correction using a GPS‐Based iterative tropospheric decomposition model. Remote Sensing of Environment, 204, 109–121. https://doi.org/10.1016/j.rse.2017.10.038

Yu, C., Li, Z., Penna, N. T., & Crippa, P. (2018). Generic atmospheric correction model for interferometric synthetic aperture radar observations. Journal of Geophysical Research: Solid Earth, 123(10), 9202–9222. https://doi.org/10.1029/2017JB015305

Yu, C., Penna, N. T., & Li, Z. (2017). Generation of real‐time mode high‐resolution water vapor fields from GPS observations. Journal of Geophysical Research: Atmospheres, 122(3), 2008–2025. https://doi.org/10.1002/2016JD025753

Zhang, Y., Gong, H., Gu, Z., Wang, R., Li, X., & Zhao, W. (2014). Characterization of land subsidence induced by groundwater withdrawals in the plain of beijing city, China. Hydrogeology Journal, 22(2), 397–409. https://doi.org/10.1007/s10040‐013‐1069‐x

Zheng, Y., & Simons, M. (2023). Ups and downs of Beverly hills, California‐analysis of surface deformation derived from InSAR: 1992 to 2023. AGU Fall Meeting Abstracts, 2023, NS32A–03.

Zhu, Y., Dortch, J. M., Massey, M. A., Haneberg, W. C., & Curl, D. (2021). An intelligent swath tool to characterize complex topographic features: Theory and application in the Teton range, licking river, and Olympus mons. Geomorphology, 387, 107–778. https://doi.org/10.1016/ j.geomorph.2021.107778