An InSAR‐GNSS Velocity Field for Iran

We present average ground‐surface velocities and strain rates for the 1.7 million km2 area of Iran, from the joint inversion of InSAR‐derived displacements and GNSS data. We generate interferograms from 7 years of Sentinel‐1 radar acquisitions, correct for tropospheric noise using the GACOS system, estimate average velocities using LiCSBAS time‐series analysis, tie this into a Eurasia‐fixed reference frame, and perform a decomposition to estimate East and Vertical velocities at 500 m spacing. Our InSAR‐GNSS velocity fields reveal predominantly diffuse crustal deformation, with localized interseismic strain accumulation along the North Tabriz, Main Kopet Dagh, Main Recent, Sharoud, and Doruneh faults. We observe signals associated with recent groundwater subsidence, co‐ and postseismic deformation, active salt diaprism, and sediment motion. We derive high‐resolution strain rate estimates on a country‐ and fault‐scale, and discuss the difficulties of mapping diffuse strain rates in areas with abundant non‐tectonic and anthropogenic signals.


Introduction
Part of the larger Alpine-Himalayan orogenic belt, convergence between the Arabian and Eurasian plates is driving active deformation and seismicity throughout Iran, concentrated in the Zagros Mountains, the Alborz, the Kopet Dagh, and the Makran subduction zone (Stern et al., 2021).Accurate geodetic estimates of ground-surface velocities and strain rates are critical to our understanding of both the localized seismic hazard, and the distribution and mechanics of deformation throughout the country.
Interferometric Synthetic Aperture Radar (InSAR) time-series analysis is a well-established technique for measuring ground deformation at the level of millimeters per year (Morishita et al., 2020), with spatial resolutions in the tens of meters, providing spatially continuous observations over hundreds of kilometres that are useful for regional continental tectonic studies (Elliott et al., 2016).Iran has been the focus of numerous InSAR-based studies over the past 20 years, capturing deformation from groundwater extraction (Babaee et al., 2020;Ghorbani et al., 2022;Haghshenas Haghighi & Motagh, 2021;Motagh et al., 2017), active salt diaprism (Barnhart & Lohman, 2012;Roosta et al., 2019), and co-and postseismic deformation (Liu et al., 2021;Penney et al., 2015;Plattner et al., 2022).Regarding interseismic strain accumulation, a common target of InSAR analysis due to the small magnitude of the signal and the implications for seismic hazard, previous Iran InSAR studies have focused on individual fault structures such as the Main Kopet Dagh Fault (Dodds et al., 2022;Walters et al., 2013), the Main Recent Fault (Watson et al., 2022), the North Tabriz Fault (Aghajany et al., 2017;Karimzadeh et al., 2013;Rizza et al., 2013), and the Doruneh Fault (Pezzo et al., 2012).Outside of Iran, since 2020 there has been a move toward unified large-scale InSAR velocity fields (e.g., Crosetto et al., 2020;Hamling et al., 2022;Ou et al., 2022;Weiss et al., 2020;Xu et al., 2021), driven partially by the abundance of open data from the European Space Agency's Sentinel-1 SAR satellite constellation.While previous country-scale velocity and strain rate fields have been derived for Iran using Global Navigation Satellite System (GNSS) station velocities (Khorrami et al., 2019;Masson et al., 2007;Vernant et al., 2004), the relatively poor spatial-resolution of Iran's GNSS network, with 10-100's km between stations, limits the ability of these studies to constrain short-wavelength tectonic and nontectonic ground motion and deformation.
We combine seven years of six and 12 day acquisitions of Sentinel-1 radar imagery with GNSS station velocities to derive a spatially near-continuous, 500 m spatial-resolution, InSAR-GNSS velocity field for a 1.7 million km 2 area encompassing the political borders of Iran.We investigate ground motion associated with interseismic strain accumulation across faults, groundwater subsidence, salt diaprism, and sediment motion with the aim to provide the current authoritative velocity field for the region.We derive high-resolution estimates of average current crustal strain rate, both on a country-scale and focused on individual major active fault structures.

Methodology
We process ∼85,000 interferograms from ∼18,500 Sentinel-1 images across 9 ascending and 10 descending tracks (Tables S1-5 in Supporting Information S1), using the semi-automated COMET-LiCSAR processing system (Lazecký et al., 2020).Images are processed in predefined frames with on average 13 bursts across each of the three Terrain Observation with Progressive Scans (TOPS) mode subswaths, with a one-burst overlap between along-track adjacent frames to aid mosaicing.On average, interferograms are formed between each image and the next three temporally adjacent images (typically 6, 12, and 18 days temporal baselines), with long-baseline interferograms generated when required to bridge data gaps or drops in coherence (i.e., snow-covered high terrain in winter).Small-baseline network lengths vary between frame, with a general maximum time-span of late 2014 to mid-2022 (see Supporting Information).We remove topographic contributions to the phase return using the 1 arcsecond SRTM Digital Elevation Model (DEM, Farr et al., 2007), and unwrap each interferogram in two dimensions using the statistical-cost network-flow algorithm (SNAPHU) version 2 (Chen & Zebker, 2000, 2001, 2002).Interferograms are multilooked by 20 in range and 4 in azimuth (∼50 m) during the processing, and are further downsampled to 0.005°× 0.005°in the subsequent time-series analysis.We mitigate tropospheric phase delay (Jolivet et al., 2014;Zebker et al., 1997) in each interferogram using the Generic Atmospheric Correction Online Service for InSAR (GACOS) (Yu et al., 2018), which performs favorably compared to other tropospheric correction techniques, but with decreasing performance for short-wavelength turbulent noise below 75 km (Murray et al., 2019).For 81/105 frames, ≥80% of interferograms show a reduction in the standard deviation of the phase of all points following the application of the GACOS correction (Figure S3 in Supporting Information S1), which we interpret as a reduction in tropospheric noise, and average line-of-sight velocities calculated from GACOS-corrected interferograms show a reduced correlation with elevation (Figure S4 in Supporting Information S1).
We generate cumulative line-of-sight displacements and average linear velocities from small-baseline networks of interferograms using LiCSBAS (Morishita et al., 2020), ensuring that there are no gaps in the displacement series for each pixel and solving the inverse problem for each pixel using least squares.We mask interferogram pixels with coherence values below 0.15, and apply manually defined masks to problematic regions.Interferograms with significant unwrapping errors, identified through triplet loop phase closure (e.g., Biggs et al., 2007), are removed based on an average loop closure threshold of 1.5 rad.Velocity uncertainties are estimated using 100 iterations of a percentile bootstrap method (Efron & Tibshirani, 1986), which results in artificially low uncertainties adjacent to the reference point (Ou et al., 2022).We calculate an experimental semivariogram for each uncertainty map and find the best-fit exponential model.The uncertainties are scaled by the ratio between the sill and the exponential model value at the distance between each pixel and the reference point (e.g., Figure S5 in Supporting Information S1).
We merge adjacent frames only along-track by removing the mode difference of the velocities (to the nearest 0.1 mm/yr) in the overlapping bursts, and then take the weighted mean of all overlapping pixels, where the weight is the inverse of the squared bootstrapped uncertainties, which improves the spatial consistency of later corrections (e.g., Figure S6 in Supporting Information S1).Grouped by track, the standard deviations of the overlap velocity differences vary from 0.6 to 3.0 mm/yr, with a mean standard deviation of 1.5 and 1.3 mm/yr for ascending and descending tracks, respectively (Figure S7 in Supporting Information S1).For comparison, Weiss et al. ( 2020) calculate along-track differences post-reference frame transformation with a standard deviation of 3.23 and 3.20 mm/yr for ascending and descending, respectively.Ou et al. (2022) minimize along-track differences using overlapping pixels, 3-D GNSS velocities, and a first-order ramp, achieving residual standard deviations of 0.3-2.1 mm/yr, with an average of 0.87 mm/yr.We achieve more consistent along-track overlaps than Weiss et al. (2020), and comparable along-track overlap differences to Ou et al. (2022).The LOS track velocities contain curved ramps in the range direction caused by the translational motion of the plate in the satellite reference frame (Stephenson et al., 2022).We mitigate these ramps by subtracting Eurasia plate velocities in the ITRF2014 plate motion model (Altamimi et al., 2016), projected into the LOS for each track, from the LOS velocities (Figures S8 and S9 in Supporting Information S1).
Average LOS velocities from LiCSBAS are initially relative to an internal reference point chosen based on the lowest root-mean-square of all loop phase closures.We transform the LOS velocities for each merged track into a Eurasia-fixed reference frame (Weiss et al., 2020;Xu et al., 2021) defined by compiled 2-D GNSS velocities from Khorrami et al. (2019).We remove 44 GNSS stations within 10 km of subsiding basins to minimize bias by these short-wavelength signals (Figure S10 in Supporting Information S1), with minimal impact on the spatial coverage of the remaining 355 stations (Figure S11 in Supporting Information S1).We interpolate the GNSS station velocities following the method of Wang and Wright (2012), employing a 0.25°triangular mesh and a smoothing factor of 10 1.5 to suppress noise in the velocities (Supporting Information, Figure S12 in Supporting Information S1).For each track, we project the East and North interpolated GNSS velocities into the satellite LOS and calculate the residual between the GNSS and InSAR velocities.We mitigate short-wavelength signals by applying a 100 km median filter, and then subtract the filtered residual velocities from the InSAR LOS velocities to transform the latter into a Eurasia-fixed reference frame (e.g., Figure S13 in Supporting Information S1).The use of a 100 km median filter limits the distortion of tectonic signals below this wavelength (e.g., interseismic strain accumulation along active faults) while mitigating spatially complex long-wavelength noise (e.g., unmitigated tropospheric, ionospheric, solid earth tides) along the full length of the track (>1,000 km) without any assumption as to the spatial structure (Supporting Information, Figures S14 and S15 in Supporting Information S1).
To aid in our interpretation of the observed ground motion, we decompose the InSAR LOS velocities into horizontal and vertical components (Ou et al., 2022;Shen, 2020).For pixels with at least one ascending (V asc ) and one descending (V dsc ) velocity, we first decompose the LOS velocities into East (V E ) and joint vertical-North (V UN ) velocity components: where θ is the local incidence angle and α is the satellite heading.V UN is then separately decomposed into V U and V N using the interpolated North GNSS velocities to constrain V N , given the relative insensitivity of Sentinel-1 InSAR to North-South motion, and the lack of vertical GNSS velocities for Iran: Given that V U is dependent upon the local incidence and azimuth angles of each given track, a maximum of four possible values of V U can be calculated if there are two ascending and two descending tracks.We calculate the weighted mean of all possible realizations of V U , using a weighting of one/σ(V U ):

Iranian Velocity Field
Using the 19 Sentinel-1 InSAR LOS track velocity fields shown in Figure 1, we generate a near-continuous East (Figure 2) and vertical (Figure 3) velocity field for a 1.7 million km 2 area encompassing almost the entirety of Iran at a pixel spacing of ∼500 m.
The decomposed East velocities show a country-scale trend of negative (westward) velocities in the West, and positive (eastward) velocities in the East, representative of the varying angle of plate convergence across Iran, as observed at large scales in previous GNSS studies (Khorrami et al., 2019).Between 45°-48°E, we observe an East velocity change between the Lesser Caucasus and the Turkish-Iranian plateau, consistent with 8-10 mm/yr of right-lateral slip along the North Tabriz Fault (NTF) and adjacent structures (Aghajany et al., 2017; Karimzadeh  (Dziewonski et al., 1981;Ekström et al., 2012) are concentrated in the major mountain belts, highlighted by shaded SRTM30 topography (Farr et al., 2007).Eurasia-fixed GNSS velocities from Khorrami et al. (2019) show the convergence of the Arabia and Eurasia plates, the boundary of which is demarked by the red line (Bird, 2003), with decreasing velocities from SW-NE (d).et al., 2013;Rizza et al., 2013;Su et al., 2016).A similar East velocity change is observed at ∼55°E, with overall positive velocities in Turkmenistan and northeastern Iran changing to negative velocities East of the South Caspian, consistent with 6-10 mm/yr of right-lateral motion along the Main Kopet Dagh Fault (MKDF) and 4-6 mm/yr left-lateral motion along the Sharoud Fault Zone (SFZ) (Dodds, 2021;Walker et al., 2021).Between 57°-60°E, we observe a small 1-2 mm/yr East velocity change across the Doruneh Fault (DF), consistent with left-lateral slip along the fault (Farbod et al., 2016;Pezzo et al., 2012;Walpersdorf et al., 2014), which is masked by the Dasht-e Kavir to the west and becomes negligible East of 60°E.We observe the same small velocity change across the Main Recent Fault (MRF) as previously detected in Watson et al. (2022), although the velocities in the western Zagros appear particularly impacted by poor coherence and phase bias.
A significant proportion of Iran's crustal deformation occurs in a North-South direction (e.g., strike-slip motion along the Minab Zendan-Palami (MZP) fault (Peyret et al., 2009) , 2019).We observe no measurable country-scale variations in the vertical velocities at the 1-2 mm/yr level, including over the major mountain ranges.However, long-wavelength vertical trends may be suppressed by the transformation of our relative LOS InSAR velocities into a Eurasia-fixed reference frame, as the GNSS velocities used to calculation the transformation (Khorrami et al., 2019) do not include vertical velocities.
We identify 13 signals associated with co-and postseismic deformation (Table S6 in Supporting Information S1) across both East and vertical velocity fields, where the estimated average velocity is distorted by non-linear postseismic deformation and/or a coseismic offset.Without the removal of both co-and postseismic deformation, care must be taken in interpreting velocities close to earthquake epicenters, as the observed ground motion may take years to decades to return to the pre-earthquake interseismic rate (Ingleby & Wright, 2017).
The vertical InSAR velocity for Iran contains a large number of subsidence signals associated with groundwater extraction (Haghshenas Haghighi & Motagh, 2021), with diameters of <10 km to >100 km and peak rates of >100 mm/yr.Most subsidence signals show seasonal variation in the cumulative displacement series correlated with seasonal variations in rainfall (Modarres & Sarhadi, 2009), with a significant variation in the strength of the seasonal term between basins.Some of the fastest subsidence signals show "pinching" signals in the East velocities, representative of inward motion of the basin sides, and care must be taken to isolate these from tectonic motion.We observe potential flexural uplift, present as "uplift halos" on the order of several millimeters per year around some subsiding basins (e.g., E-E' in Figure 3), caused by unloading of the lithosphere (Amos et al., 2014;Chaussard et al., 2017;Hammond et al., 2016).
In the southeastern Zagros Mountains, active salt diaprism of the Hormuz salt layer (Kent, 1979;Talbot, 1998) has previously been measured at a regional-scale using Envisat by Barnhart and Lohman (2012), who identify salt-related ground motion signals at 20 locations, and a lack of measurable ground motion at a further 30 diapirs.We confirm many of the same active versus inactive assessments in our East velocities, with several diapirs showing either motion on previously 'inactive' diapirs, or vice versa.While these changes could reflect accelerations or deceleration of deformation at these diapirs, with the small-baseline networks from Barnhart and Lohman (2012) spanning 2004-2010versus 2014-2022 in this study, they may also represent a greater signal-tonoise ratio in our velocity estimates, primarily as a result of the vastly shorter repeat time of Sentinel-1 compared to Envisat.Further work is required to map spatial and temporal variations in active salt diaprism across the Zagros, with a robust comparison to Barnhart and Lohman (2012).

Geophysical Research Letters
10.1029/2024GL108440 We observe elevated East and vertical velocities throughout the Dasht-e Kavir in the Central Iranian Plateau, broadly correlated with playa lakes, which may capture apparent ground motion from sediment flow and deposition.However, these signals may be heavily influenced by phase bias (De Zan et al., 2015;De Zan & Gomba, 2018), which positively correlates with decreasing interferogram temporal baseline (Ansari et al., 2020;Purcell et al., 2022), and which has been observed around salt lakes in arid regions (Zheng et al., 2022).Future InSAR work in Iran must account for phase bias, either through empirical corrections (Maghsoudi et al., 2022;Zheng et al., 2022) or avoidance of short-temporal-baseline interferograms (Purcell et al., 2022), even in regions previously thought to be minimally impacted by phase bias due to their high aridity.

Strain Rates
The decomposed East velocities suggest that strain is relatively diffuse across most of Iran, with a small number of faults exhibiting localized interseismic strain accumulation that is observable in Sentinel-1 InSAR time-series.
We calculate components of the strain rate tensor directly from the decomposed East velocities and the interpolated North GNSS velocities (Supporting Information), using a sliding circular median filter to suppress random noise and short-wavelength non-tectonic signals in the East velocities (Ou et al., 2022), and potential mesh artifacts in the North velocities.The optimal window size for isolating interseismic strain accumulation along active faults is a product of the signal-to-noise ratio (Ou et al., 2022), the spatial wavelength of the deformation signal, and the relative adjacency of signals with similar magnitudes.We experiment with a range of filter window widths (4-, 10-, 20-, 40-, 60-, 80-, and 100 km) and apply a manually defined mask for co-and postseismic signals and subsidence to the East and North velocities, opting not to remove the ground motion signals in the Dasht-e Kavir so as to highlight their impact on attempts to isolate the underlying low-magnitude tectonic strain rate (Masson et al., 2014).
At the 80 km scale, we observe elevated strain rates (≥50 nst/yr) throughout the Zagros Mountains, at the western edge of the Makran, around the edges of the Dasht-e Lut, and along the NTF, MKDF, and SFZ (Figure 4).In the northwestern Zagros (∼47°E), we observe potential along-strike postseismic strain from the Sarpol-e Zahab earthquake, while in the southeastern Zagros (52°-56°E), diffuse strain rates may reflect the impact of the Hormuz Salt forming a detachment horizon.The Central Iranian Plateau exhibits a lower overall strain rate, matching the expected lack of internal deformation based on seismicity (Karasözen et al., 2019;Raeesi et al., 2017), with the exception of the Dasht-e Kavir where potential non-tectonic surface motion masks any underlying tectonic signal.We observe low strain rates in Turkmenistan, where the Turan platform is commonly assumed to be part of the stable continental interior (e.g., Jackson et al., 2002) and along the western edge of Afghanistan, The central Makran (58°-63°E) exhibits a low average strain rate, potentially indicating a lower degree of coupling on the subduction interface, although this is in disagreement with previous studies (Lin et al., 2015;Penney et al., 2017) and may be explained by poor North-South sensitivity in the InSAR and a suppression of long-wavelength vertical signals by the referencing transformation.In contrast, the western Makran shows rapid and diffuse strain accumulation related to the Oman indenter and N-S motion along the Minab Zendan-Palami fault zone.
A 60 km median filter window width (as proposed by Ou (2020)) proves effective at mitigating shortwavelength noise and isolating the fault signal across the MKDF and SFZ.For filter widths between 40 and 80 km, the MKDF exhibits elevated strain rates of 100-150 nst/yr.However, the masked area to the North of the MKDF associated with subsidence in the Karakum Canal Irrigated Zone (Dodds, 2021), potentially causes a shift of the strain rate signal to the North of the MKDF at higher filter widths, with a 40 km filter (profile C-C') keeping the signal centered on the mapped fault trace.Irrespective of this shift, the elevated strain rate signal along the MKDF follows the mapped fault trace between 56°-57°E, but shifts 30-40 km southwest of the fault at the northwestern end (55°-56°E), potentially onto a sub-parallel splay fault (Dodds, 2021).The SFZ exhibits higher strain rates, with a range of 150-200 nst/yr across the fault.At smaller filter window sizes (10 km), we can identify potential strain localisation on two adjacent faults, information that is lost at higher filter window sizes (20-40 km).
We observe elevated (300-400 nst/yr at a 4 km length-scale, down to 100-250 nst/yr at a 20 km length-scale) strain rates associated with the mapped trace of the Doruneh Fault, with a short-spatial-wavelength indicative of either shallow interseismic strain accumulation, or near-surface creep, although the latter cannot be determined given the 500 m pixel spacing of the velocity field.The short-wavelength of the strain signal requires a small 4 km filter window to avoid suppressing the tectonic signal, resulting in significant off-fault noise of a comparable magnitude.The NTF exhibits a broad (∼100 km) zone of elevated strain along the length of the fault, although much of the signal is lost to both incoherence and masking of the 2019 Bozghush earthquake signal (Isik et al., 2021;Yang et al., 2020).
When interpreting strain rates calculated from filtered InSAR velocities, care must be taken as the use of a median filter produces Gibbs-ringing artifacts (Kellner et al., 2016) -oscillating sub-parallel signals adjacent to the main tectonic signals -which may be over interpreted as distributed strain across multiple faults.Additionally, the choice of filter window size determines the separation distance at which adjacent tectonic fault signals can be isolated from one-another.

Conclusion
We have produced the first country-scale combined InSAR-GNSS velocity for Iran, encompassing an area of 1.7 million km 2 at a ∼500 m pixel spacing, and derived high-resolution estimates of average current crustal strain rate.Our East velocities highlight a complex mix of tectonic and non-tectonic signals, with localized interseismic strain accumulation along the North Tabriz, the Main Kopet Dagh, the Sharoud, and the Doruneh faults.While the Central Iranian Plateau exhibits low internal deformation, significant regions are masked by potential sediment motion and phase bias throughout the Dasht-e Kavir, a region otherwise assumed to be of low risk from phase bias due to the high aridity.The vertical velocity field highlights an abundance of subsidence signals associated with groundwater extraction, with associated horizontal motion and potential flexural uplift halos for the largest of these signals.High-resolution North-South velocity estimates, both from burst overlaps and the future SAR satellite missions, will be critical for measuring the complex and diffuse deformation observed throughout Iran.Figures were made using the Generic Mapping Tools (Wessel et al., 2013).We thank Qi Ou for strain rate discussions.Andrew Watson is supported through a PhD studentship from the Royal Society (RGF\R1\180076).This work is supported by the UK Natural Environment Research Council (NERC) through the Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET, http://comet.nerc.ac.uk).This work is also supported by NERC through the Looking into the Continents from Space (LiCS) large Grant (NE/K010867/1).John Elliott gratefully acknowledges support from the Royal Society through a University Research Fellowship (UF150282 and URF\R\211006).

Figure 1 .
Figure 1.Sentinel-1 line-of-sight (LOS) InSAR-derived average surface velocities for COMET-LiCSAR-defined frames (dashed black polygons) in ascending (a) and descending (b) look directions, where positive values (red) indicate apparent motion toward the satellite.Seismotectonic overview of Iran (c), with the major tectonic units labeled and fault traces from the Global Earthquake Model (GEM) active fault database(Styron & Pagani, 2020) shown in gray.Focal mechanisms from the Global Centroid Moment Tensor (GCMT) catalog(Dziewonski et al., 1981;Ekström et al., 2012) are concentrated in the major mountain belts, highlighted by shaded SRTM30 topography(Farr et al., 2007).Eurasia-fixed GNSS velocities fromKhorrami et al. (2019) show the convergence of the Arabia and Eurasia plates, the boundary of which is demarked by the red line(Bird, 2003), with decreasing velocities from SW-NE (d).

Figure 3 .
Figure 3. Decomposed Vertical InSAR-GNSS velocities, with positive (red) values indicating upward motion (see Figure 2 caption for repeated elements).Example cumulative displacement series are given for four subsiding basins (a -Asadabad, (b) Isfahan, (c) Damghan, (d) Rafsanjan), with the average linear velocity (V LOS ) estimated by LiCSBAS given in the top left.Profile E-E' contains potential flexural uplift (shaded regions in inset) either side of a subsiding basin.Fault traces from the GEM active fault database are shown in gray.

Figure 4 .
Figure 4. Second invariant of the strain rate tensor derived from the decomposed East velocities and the interpolated GNSS North velocities for appropriate median filter widths, with fault traces from the GEM active fault database shown in black and gray.(a) Strain rate for all of Iran (80 km).(b) Long-wavelength elevated strain over the NTF (100 km).(c) Elevated strain across the MKDF and SFZ (60 km).The dashed black line defines a potential splay fault implied by the strain rates.(d) Shortwavelength elevated strain across the DF (4 km).For six across-fault profiles (A-F), we show average profiled strain rates (1 km long bins, 5 km wide swath) for different filter window widths.Karakum Canal Irrigated Zone (KCIZ), Quchan-Bakhardan Fault Zone (QBFZ).
Li et al., 2021)al., 2014)faults(Walpersdorf et al., 2014)), to which we have very limited sensitivity in conventional line-of-sight Sentinel-1 InSAR observations.Future work should incorporate measurements of North-South ground motion from Sentinel-1 Burst-Overlap Interferometry (e.g.,Li et al., 2021), and later from the ESA Harmony Satellite mission (López-Dekker et al.