# Large-Scale Interseismic Strain Mapping of the NE Tibetan Plateau From Sentinel-1 Interferometry

Q. Ou, S. Daout, J. Weiss, L. Shen, M. Lazecký, T. Wright, B. Parsons

► To cite this version:

Q. Ou, S. Daout, J. Weiss, L. Shen, M. Lazecký, et al., Large-Scale Interseismic Strain Mapping of the NE Tibetan Plateau From Sentinel-1 Interferometry. Journal of Geophysical Research: Solid Earth, 2022, 127 (6), pp.1-29. 10.1029/2022JB024176. [https://doi.org/10.1029/2022JB024174]

HAL Id: hal-04494144

https://hal.science/hal-04494144v1

Submitted on 10 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2022JB024176

**Special Section:**

100-year Anniversary of the Great 1920 Haiyuan Earthquake: What have We Learnt on Large Continental Earthquakes and Faults?

**Key Points:**

* We present new methods to mosaic line-of-sight velocity frames, derive Cartesian velocities, and produce high-resolution strain-rate maps

* We construct 440,000 km² maps of east-west, vertical velocities, and strain rates over NE Tibet from Sentinel-1 interferometry

* We quantify strain partitioning, measure creep rate, identify unmapped structures and climatic, hydrological and anthropogenic signals

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

Q. Ou,
qi.ou@earth.ox.ac.uk

**Citation:**

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, e2022JB024176. https://doi.org/10.1029/2022JB024176

Received 11 FEB 2022
Accepted 24 MAY 2022

**Author Contributions:**

**Conceptualization:** T. J. Wright, B. E. Parsons

**Data curation:** Q. Ou, M. Lazecký

**Formal analysis:** Q. Ou

**Funding acquisition:** T. J. Wright, B. E. Parsons

© 2022. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

## Large-Scale Interseismic Strain Mapping of the NE Tibetan Plateau From Sentinel-1 Interferometry

Q. Ou¹ ID, S. Daout¹,² ID, J. R. Weiss³,⁴,⁵ ID, L. Shen³ ID, M. Lazecký³ ID, T. J. Wright³ ID, and B. E. Parsons¹ ID

¹COMET, Department of Earth Sciences, University of Oxford, Oxford, UK, ²CRPG, Université de Lorraine - CNRS, UMR 7358, rue N. D. des Pauvres, Vandœuvre les Nancy, France, ³COMET, School of Earth and Environment, University of Leeds, Leeds, UK, ⁴Institute of Geosciences, University of Potsdam, Potsdam, Germany, ⁵Present at NOAA/NWS Pacific Tsunami Warning Center, Honolulu, HI, USA

**Abstract** The launches of the Sentinel-1 synthetic aperture radar satellites in 2014 and 2016 started a new era of high-resolution velocity and strain rate mapping for the continents. However, multiple challenges exist in tying independently processed velocity data sets to a common reference frame and producing high-resolution strain rate fields. We analyze Sentinel-1 data acquired between 2014 and 2019 over the northeast Tibetan Plateau, and develop new methods to derive east and vertical velocities with ~100 m resolution and ~1 mm/yr accuracy across an area of 440,000 km². By implementing a new method of combining horizontal gradients of filtered east and interpolated north velocities, we derive the first ~1 km resolution strain rate field for this tectonically active region. The strain rate fields show concentrated shear strain along the Haiyuan and East Kunlun Faults, and local contractional strain on fault junctions, within the Qilianshan thrusts, and around the Longyangxia Reservoir. The Laohushan-Jingtai creeping section of the Haiyuan Fault is highlighted in our data set by extremely rapid strain rates. Strain across unknown portions of the Haiyuan Fault system, including shear on the eastern extension of the Dabanshan Fault and contraction at the western flank of the Quwushan, highlight unmapped tectonic structures. In addition to the uplift across most of the lowlands, the vertical velocities also contain climatic, hydrological or anthropogenic-related deformation signals. We demonstrate the enhanced view of large-scale active tectonic processes provided by high-resolution velocities and strain rates derived from Sentinel-1 data and highlight associated wide-ranging research applications.

**Plain Language Summary** The new-generation radar satellites permit the derivation of crustal velocities in fine resolution over continental scales. With the technical advancements developed in this study, we reveal unseen details of the crustal deformation over the earthquake-prone Hexi Corridor in the northeast corner of the Tibetan Plateau. In particular, we observe concentrated elastic loading along, and at branch points of, major earthquake-generating faults. The ability to monitor such phenomena is important for understanding future seismic hazard as such geometrical complexities along faults are often associated with earthquake triggering and termination. We also draw attention to previously unknown structures that are rapidly deforming. This new and enhanced view of the local tectonic setting could improve local seismic hazard assessment. Last but not least, we observe climatic, hydrological and anthropogenic signals in the vertical velocity field, which are related to permafrost thawing, blocked river drainage, mining, damming and the extraction of groundwater for farming. Overall, we demonstrate the wealth of information that can be derived from the rapidly growing space-born Earth observation data, which are destined to play an important role in shaping a more resilient world for the future.

## 1. Introduction

Over recent decades, Interferometric Synthetic Aperture Radar (InSAR) has emerged as a powerful geodetic tool for imaging crustal deformation on the continents (e.g., Burgmann et al., 2000; Elliott et al., 2016). Numerous studies have demonstrated InSAR's ability to reveal the presence of active faults that were not previously recognized (e.g., Daout et al., 2019; Wicks et al., 2013), highlight strain concentration (Weiss et al., 2020), identify creeping sections of faults (Cavalié et al., 2008; Jolivet et al., 2012; Rousset et al., 2016), and constrain slip parameters and frictional properties of faults (Jolivet et al., 2013, 2015; Zhou et al., 2018). This information is fundamental to seismic hazard assessment and understanding crustal dynamics.
**Investigation:** Q. Ou

**Methodology:** Q. Ou, S. Daout, J. R. Weiss, L. Shen

**Project Administration:** T. J. Wright, B. E. Parsons

**Software:** S. Daout, J. R. Weiss

**Supervision:** S. Daout, T. J. Wright, B. E. Parsons

**Validation:** S. Daout

**Visualization:** Q. Ou

**Writing – original draft:** Q. Ou

**Writing – review & editing:** Q. Ou, S. Daout, J. R. Weiss, L. Shen, M. Lazecký, T. J. Wright, B. E. Parsons

The global coverage and the regular and frequent acquisition of SAR measurements provided by the Sentinel-1A and Sentinel-1B satellites, launched by the European Space Agency in April 2014 and April 2016 respectively, opened the prospect of mapping high-resolution crustal velocities and strain rates worldwide. Various processing systems have been developed to produce interferograms and time series from the continuously growing SAR data sets, such as NASA Jet Propulsion Laboratory's Advanced Rapid Imaging and Analysis (ARIA) (Bekaert et al., 2020), COMET's Looking Into Continents from Space with Synthetic Aperture Radar (LiCSAR) and LiCSBAS suites (Lazecký et al., 2020; Morishita et al., 2020) and ForM@Ter Solid Earth data and services center's ForM@Ter LArge-scale multi-Temporal Sentinel-1 InterferoMetry (FLATSIM) (Thollard et al., 2021).

For automatic large-scale high-resolution velocity mapping projects, subdividing the target area into smaller patches for batch processing is an efficient way of using finite computing power, memory, and disk space. In contrast to processing InSAR data in long tracks (e.g., Thollard et al., 2021), a frame-based processing approach offers the flexibility of gradually growing the mapping area to eventually cover the full continent. It also allows unwrapping challenges to be tackled and interferogram networks to be refined on a frame by frame basis (e.g., Weiss et al., 2020). In the case of LiCSAR, acquisitions are organized into ~250 × 250 km² frame units. Yet, multiple challenges exist in combining the independently processed velocity frames into large-scale velocity fields and producing high-resolution strain rate fields from the velocities. Mosaiced line-of-sight (LOS) velocities can have large mismatches between frames along track (Weiss et al., 2020), which are only zero if identical networks of interferograms are used for all frames. The high-resolution velocity data often need to be downsampled for Cartesian velocity decomposition and strain rate calculation due to the associated computational cost (H. Wang & Wright, 2012; Xu et al., 2021). The resultant strain rate fields can appear overly smooth (Song et al., 2019; Weiss et al., 2020). This study aims to tackle these technical challenges in order to derive large-scale high-resolution velocity and strain rate results without compromising the high quality of the Sentinel-1 data.

We investigate the rapidly deforming NE Tibetan Plateau (Figure 1a), at the leading edge of the northward expansion of the Tibetan Plateau driven by the ongoing collision of India with Eurasia (England & Houseman, 1985; England & Molnar, 1997; Flesch et al., 2001; Molnar & Tapponnier, 1975; Pichon et al., 1992; Tapponnier et al., 2001; Yuan et al., 2013). This area hosted 20 $M_W$ > 6.5 earthquakes in the past century, including the $M_W$ 7.9 1920 Haiyuan Earthquake (Deng et al., 1986; IGCEA & NBCEA, 1990; Ou et al., 2020; Ren et al., 2016; Xu et al., 2019). Extensive research has been carried out to study the kinematics (P. Zhang et al., 1988; W. Zheng et al., 2013), mechanics (Deng et al., 1984; Gaudemer et al., 1995), and seismic hazard (Liu-zeng et al., 2015; Xiong et al., 2010) of the major faults in the region, such as the Haiyuan, Kunlun and West Qinling Faults and the Qilianshan thrusts. However, how strain, the prerequisite for earthquakes, partitions between the major and minor faults, and off-fault areas remains unclear. Geological slip rates have mostly been estimated for the fast-slipping major faults (Kirby et al., 2007; C. Li et al., 2009; J. Li et al., 2016; Yao et al., 2019; Shao et al., 2020). The GNSS network in the area is too sparse to resolve the degree of strain localisation (Gan et al., 2007; Liang et al., 2013; G. Zheng et al., 2017; M. Wang & Shen, 2020). Previous InSAR studies in the area had limited spatial coverage (Cavalié et al., 2008; Daout et al., 2016; Jolivet et al., 2012; Qiao et al., 2021; Song et al., 2019).

We aim to provide the first high-resolution regional overview of present-day deformation of the NE Tibetan Plateau. The large-scale InSAR velocity and strain rate fields resulting from this study provide the fundamental data sets for studying strain partitioning, crustal rheology, and seismic hazard, for testing competing theories about how the continental lithosphere deforms, and for monitoring climatic, hydrological, and anthropogenic processes.

In the following sections, we first describe the processing workflow for generating LOS velocity and uncertainty maps for each frame (Section 2). Then, we introduce new methods for mosaicing and georeferencing the frame maps using GNSS data (Section 3), decomposing the LOS velocity tracks into regional east ($V_E$) and vertical ($V_U$) velocity maps (Section 4) and producing low-noise, high-resolution strain rate maps (Section 5). Finally, we highlight the main features observed in the resultant $V_E$, $V_U$ and strain rate maps, and discuss the implications for regional tectonics, fault mechanics and seismic hazard (Section 6).
<figure><img src="image_3.png" /><figcaption>Figure 1.</figcaption></figure>
## 2. InSAR Data Processing and Time Series Analysis

### 2.1. Data

Sentinel-1 images over the NE Tibetan plateau were acquired on ascending and descending passes every 24 days from late 2014, and every 12 days from January 2017. We processed 10 ascending frames and 13 descending frames (as defined by LiCSAR, Figures 1b and 1c) of Sentinel-1A and 1B level-1 Single Look Complex (SLC) products acquired between October 2014 and December 2019. The frames have 65–110 (average 95) acquisition epochs over time spans of 3.2–5.2 (average 4.6) years (Figure S1–S2 in Supporting Information S1). Based on statistical reduction of noise, the uncertainty of the LOS velocities derived from time series data, $\sigma_v$, is expected to be given by:

$$ \sigma_v \approx \frac{2\sqrt{3}\sigma_e}{\sqrt{NT}} \qquad (1) $$

where $\sigma_e$ is the uncertainty of cumulative displacement per epoch, N is the number of epochs and T is the time span in years (J. Zhang et al., 1997; Morishita et al., 2020). Assuming $\sigma_e \sim 10$ mm as in Morishita et al. (2020) obtained from a comparison of InSAR and GNSS time series, the number of epochs and lengths of the time span of data we process would theoretically give an approximate lower bound on uncertainty, $\sigma_v$, of 0.77 mm/yr in the final velocity maps.

### 2.2. Interferogram Generation

Standard 4 × 20 multilooked interferograms of ~100 m resolution were generated from the Sentinel-1 imagery with the LiCSAR automatic processing package (Lazecký et al., 2020). By default, LiCSAR generates interferograms between each epoch and three consecutive epochs both forwards and backwards in time. We changed the limit to six consecutive epochs to create greater redundancy, and added interferograms with 6-month and 9-month temporal baselines to provide further connections and hence reduce the impact of the potential biases from the short temporal baseline interferograms (Ansari et al., 2020; Daout et al., 2020). This gave an average of 540 interferograms for each frame (total of 12,480). Interferograms were filtered and unwrapped using version two of the statistical-cost, network-flow phase-unwrapping (SNAPHU) algorithm (C. W. Chen & Zebker, 2002), as described in Lazecký et al. (2020). The 23 frames have on average 86% of pixels unwrapped (Figure S1 and S2 in Supporting Information S1), as expected from the high coherence of the data aided by the stable orbital control of the Sentinel satellites and short temporal baselines. Subsequent processing was carried out with the interferograms in radar coordinates.

### 2.3. Atmospheric Correction

To correct for tropospheric delays in the unwrapped interferograms, we used the Generic Atmospheric Correction Online Service (GACOS) model (Yu et al., 2018), which performs similarly if not better than the reanalysis model from European Center for Medium-Range Weather Forecasts (ERA5) in the Tibetan Plateau and in most regions of the world (Y. Wang et al., 2021; Yu et al., 2021). The GACOS models per epoch were converted to radar coordinates, back-projected to the line-of-sight (LOS) direction, and differenced in the same order as the interferometric pairs, before being subtracted from the unwrapped interferograms (Figures 2a, 2b and 2d).

The reference area was chosen manually for each frame in an area that is expected to be temporally stable and away from active faults, vegetation, agriculture and desert. The reference window size was set to be 400 × 400 pixels², which is ~30 × 30 km². This size is small enough to avoid tectonic signals and large enough to average over high-frequency turbulent atmospheric noise. The average pixel value in the reference window is subtracted from the interferograms and their corresponding GACOS models, so that both data sets have the same local reference.

**Figure 1.** (a) Seismo-tectonic map of the NE Tibetan Plateau with seismicity from ISC-GEM (ISC, 2019) and GCMT (Dziewonski et al., 1981; Ekström et al., 2012) catalogs, GPS velocity vectors from M. Wang and Shen (2020), topography from the 30-arcsecond SRTM digital elevation model (DEM) (Farr et al., 2007) and faults adapted from Yuan et al. (2013) W. Zheng et al. (2016) and (b) Li et al. (2020). Red and blue polygons are the outlines of areas covered by ascending and descending Sentinel-1 Interferometric Synthetic Aperture Radar frames, respectively. (b and c) show the locations of ascending and descending frames respectively, with COMET-LICSAR frame names labeled.
<figure><img src="image_4.png" /><figcaption>Figure 2.</figcaption></figure>
We evaluated the effectiveness of the GACOS correction using both the best-fit slope and Pearson correlation coefficient of the scatter plots between GACOS and interferogram phases (Figure 2c). The 23 frames have correlations between 0.50 and 0.78 and slopes between 0.37 and 0.71 (Figure S3 and S4 in Supporting Information S1), suggesting that GACOS is better at predicting the spatial pattern of the tropospheric effect than predicting its amplitude (L. Shen et al., 2019). Interferograms with larger amplitude ranges tend to have slopes and correlations closer to 1 (Figure S3 and S4 in Supporting Information S1), suggesting the most significant atmospheric delays in the large-amplitude interferograms have been the most effectively corrected.

## 2.4. Ramp Removal

A planar ramp was then removed from each interferogram by solving for a correction term for each epoch using a network approach following Biggs et al. (2007) and Maubant et al. (2020) (Figures 2e–2g, Text S1 in Supporting Information S1). This step reduces any remaining long-wavelength signals that might be associated with unmodelled tropospheric delays, ionospheric phases (Gomba et al., 2017), orbital inaccuracies (Fattahi & Amelung, 2014), solid earth tide (Xu & Sandwell, 2020) and any long wavelength tectonics. Flattening all interferograms at this step also helps reveal unwrapping error during network refinement.

## 2.5. Network Refinement

Time series analysis was performed using the NSBAS program developed by López-Quiroz et al. (2009) and Doin et al. (2015). We used the misfit maps per interferogram (Figure 2h), per epoch (Figure 2i) and their root-mean-squared (RMS) values, $\phi_{RMS_{IFG}}$ and $\phi_{RMS_{epoch}}$ to identify interferograms with unwrapping error. Regions of high values in the $\phi_{RMS_{pixel}}$ map for the entire network (Figure 2j) also indicate the presence of interferograms with unwrapping error. These erroneous interferograms were discarded during the preliminary inversion, which was performed at low resolution and without temporal smoothing or automatic corrections. If the resulting network becomes poorly connected, we made additional interferograms with longer temporal baselines using LiCSAR to build up the network for another preliminary inversion. This process was repeated until all the interferograms in a well-connected network (with at least three interferograms linking every epoch) have $\phi_{RMS_{IFG}}$ lower than ~2 rad, and no more clear unwrapping error appear in the overall $\phi_{RMS_{pixel}}$ map.

Then, we remove all 12-day interferograms from the networks for the final time series inversion to mitigate the potential bias from interferograms with short temporal baselines (Ansari et al., 2020). As the 12-day interferograms are only present in the post-2017 time series, removing them means the entire network from 2014 to 2019 is now made of interferograms with temporal baseline of 24 days and above. Removing them only before the final inversion allows us to capitalize on the coherent pixels in those 12-day interferograms for better ramp removal and error identification in the preliminary inversions. The final networks remain connected and have an average of 430 interferograms per frame (Figure S1 and S2 in Supporting Information S1). Similar to Daout et al. (2020) and Weiss et al. (2020), we find the bias caused by the short-period (12 and 24 days) interferograms negligible when the network is well-connected and contains many long-period (3 and 6 months) interferograms.

## 2.6. Time Series Analysis

The final network was then inverted at the same resolution as the multilooked interferograms (100 m) with temporal smoothing and automatic correction (López-Quiroz et al., 2009) (Text S2). The resultant $\phi_{RMS_{pixel}}$ map

**Figure 2.** (a–g) Atmospheric delay and ramp removal for an example of an unwrapped interferogram generated by Looking Into Continents from Space with Synthetic Aperture Radar between epochs 25 December 2016 and 18 January 2017 from the frame 033D_05304_131313. The gray boxes in panels (a–g) outline the reference window manually chosen in a stable non-deforming region. The mean value in the reference windows are subtracted from both (a) the unwrapped interferogram and (b) the Generic Atmospheric Correction Online Service (GACOS) atmosphere model back-projected to the line-of-sight (LOS) direction. (c) Scatter plot between corresponding pixel values in (a) and (b) to calculate the slope of the best-fit red dashed line and the correlation of the interferogram and GACOS model for statistical analysis (Figures S3 and S4). (e) The best-fit ramp, estimated from (d) the residual after GACOS correction, can be different from (f) the final ramp subtracted from (g) to give (h) the flattened interferogram because the ramp in (f) is obtained through a network inversion to ensure phase closure around interferogram loops (h and i) Unwrapping error revealed by (h) the error map generated between input and reconstructed interferograms 21 August 2016–18 January 2017, (i) the root-mean-squared (RMS) errors averaged across all interferograms associated with epoch 21 August 2016 and (j) the RMS errors averaged across the preliminary network. All interferograms associated with epoch 21 August 2016 were found to contain unwrapping error and were removed from the network (k) Improved global error map after network refinement (l and m) LOS velocities and (o and n) uncertainties in radar and geographical coordinates, respectively, with gray box in (o) highlighting low uncertainties around the reference window.
was visually checked to confirm there are no remaining unwrapping error (Figure 2k). We inverted a linear velocity, V(LOS), from the time series using weights derived from $\phi_{RMS_{epoch}}$, deviation of each cumulative displacement map from a plane, and deviation of the time series of each pixel from a linear velocity model (Daout et al., 2017, 2018) (Figure 2l and 2o, Text S3). The uncertainty, $\sigma(LOS)$, takes into account the misfit, the model variance (Tarantola, 2005) and the numbers of degrees of freedom in the velocity inversion (difference between the numbers of data points and model variables) (Equation 2 in Supporting Information of Daout et al. (2017)) The resultant linear velocity and uncertainty maps were then cleaned with the $\phi_{RMS_{pixel}}$ map with a threshold of 0.5 radian to only retain the highest quality pixels. Figures 2m and 2n show the geocoded V(LOS) (positive toward the satellite) and $\sigma(LOS)$ maps with their units converted to mm/yr.

### 3. Frame Mosaicing

When the 23 LOS velocity frames are placed in the same regional map, frame and track boundaries are discernible in many places (Figure S5). Apart from the local referencing that caused constant offsets between frames, and the change of incidence angles that caused offsets at track overlaps, the ramps removed from the interferograms, especially from frames containing velocity steps across faults, also tilted the LOS velocity frames relative to each other.

To mosaic the frames, Weiss et al. (2020) applied a second-order polynomial to fit the InSAR LOS to the interpolated GNSS LOS, following the method of Hussain et al. (2018). The resultant offsets between frame overlaps along- and across-track had standard deviations between 3.08 and 3.72 mm/yr (Weiss et al., 2020). These values are significantly larger than the ~1 mm/yr level (0.77 × $\sqrt{2}$; see Section 2.1) expected from a 5-year time series. For frame offsets along track, the standard deviation should be even smaller given the largely shared atmospheric noise at the frame overlaps (zero if the networks are identical).

#### 3.1. Joint Inversion of InSAR Frame Overlaps and GNSS Data Points

We aim to develop a method that minimizes the LOS differences between frame overlaps along track without using interpolated GNSS velocities for frame tying. As our frame and track offsets did not show compelling evidence for the need of higher order adjustments, and considering that higher-order ramps could potentially remove tectonic signals, we invert for a planar ramp per frame to minimize alteration of the InSAR data.

Our method stems from a more advanced method developed by L. Shen (2020) to mosaic 12 ascending and descending LOS velocity tracks with limited GNSS data. We have more (23) velocities frames to mosaic, thus more (69) ramp parameters to invert for. However, we benefit from a better distribution of GNSS velocities, which allow us to mosaic frames track by track. Our joint inversion is constructed as follows to invert for a planar ramp per frame that fits the LOS differences both between overlapping InSAR frames and between InSAR and GNSS points (Figure S6 in Supporting Information S1):

$$ \begin{bmatrix} \frac{1}{\sqrt{\sigma_{11}^2 + \sigma_{12}^2}}(i1 - i2) \\ \vdots \\ \frac{1}{\sqrt{\sigma_{12}^2 + \sigma_{13}^2}}(i2 - i3) \\ \vdots \\ \frac{\sqrt{g/h}}{\sqrt{\sigma_{11}^2 + \sigma_{G1}^2}}(I1 - G1) \\ \vdots \\ \frac{\sqrt{g/h}}{\sqrt{\sigma_{12}^2 + \sigma_{G2}^2}}(I2 - G2) \\ \vdots \\ \frac{\sqrt{g/h}}{\sqrt{\sigma_{13}^2 + \sigma_{G3}^2}}(I3 - G3) \end{bmatrix} = \begin{bmatrix} \frac{1}{\sqrt{\sigma_{11}^2 + \sigma_{12}^2}}(x_{12}) & y_{12} & 1 & -x_{12} & -y_{12} & -1 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \frac{1}{\sqrt{\sigma_{12}^2 + \sigma_{13}^2}}(0) & 0 & 0 & x_{23} & y_{23} & 1 & -x_{23} & -y_{23} & -1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \frac{\sqrt{g/h}}{\sqrt{\sigma_{11}^2 + \sigma_{G1}^2}}(x_1) & y_1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \frac{\sqrt{g/h}}{\sqrt{\sigma_{12}^2 + \sigma_{G2}^2}}(0) & 0 & 0 & x_2 & y_2 & 1 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \frac{\sqrt{g/h}}{\sqrt{\sigma_{13}^2 + \sigma_{G3}^2}}(0) & 0 & 0 & 0 & 0 & 0 & x_3 & y_3 & 1 \end{bmatrix} \begin{bmatrix} a_1 \\ b_1 \\ c_1 \\ a_2 \\ b_2 \\ c_2 \\ a_3 \\ b_3 \\ c_3 \end{bmatrix} \quad (2) $$
where *l* in 1, 2 and 3 represents successive InSAR frames along track, *il* are InSAR LOS from the *l*th frame in the frame overlap with uncertainties *σ*<sub>il</sub>, *Il* are InSAR LOS in the *l*th frame at GNSS locations with uncertainties *σ*<sub>Il</sub>, *Gl* are GNSS velocities projected into the line of sight in the *l*th frame with uncertainties *σ*<sub>Gl</sub>; *a*<sub>l</sub>, *b*<sub>l</sub>, *c*<sub>l</sub> are the ramp parameters to be inverted for each frame; *x* and *y* are the coordinates of the corresponding pixels in the track, with subscripts 12 and 23 indicating overlaps between frames 1 and 2 and between frames 2 and 3, respectively, *x*<sub>i</sub> and *y*<sub>i</sub> indicate coordinates of the GNSS points in the *l*th frame; *g* is the total number of InSAR-InSAR offsets from all frame overlaps, and *h* is the total number of InSAR-GNSS overlaps from all frames.

For the top two rows representing InSAR frame overlap differences, the InSAR LOS values (*il*) are taken from the *V*(LOS) maps (Figure 2m) and the *σ*<sub>il</sub> are taken from the *σ*(LOS) maps (Figure 2n). For the bottom three rows representing InSAR-GNSS offsets, the InSAR LOS velocities (*Il*) are chosen to be the weighted average InSAR LOS velocity in a window surrounding the GNSS location, using *σ*(LOS) as weights. This is designed to handle local InSAR variability and avoid wasting GNSS points at empty InSAR pixels. We initialized the size of the window as a 5 × 5 grid, and let it grow in steps of 2 until it finds at least 10 non-empty pixels from InSAR. All stations found enough InSAR pixels within a window size of 61 × 61, which is equivalent to an area of ~4.3 × 4.3 km², smaller in dimension than that over which short wavelength tectonic motions like creep would typically occur. The uncertainty associated with this weighted mean is set to be the standard deviation of InSAR LOS pixels in the window.

The GNSS LOS velocities, *Gl*, and their associated uncertainties, *σ*<sub>Gl</sub>, are calculated as:

$$Gl = -V_E \cos(\phi) \sin(\theta) + V_N \sin(\phi) \sin(\theta) + V_U \cos(\theta) \quad (3)$$

$$\sigma_{Gl} = \sqrt{(\sigma_E \cos(\phi) \sin(\theta))^2 + (\sigma_N \sin(\phi) \sin(\theta))^2 + (\sigma_U \cos(\theta))^2} \quad (4)$$

where $V_E$, $V_N$ and $V_U$ are the GNSS velocity components in the east, north and vertical directions, and $\sigma_E$, $\sigma_N$ and $\sigma_U$ are their associated uncertainties. The horizontal GNSS velocities are in a fixed-Eurasia reference frame (M. Wang & Shen, 2020) and the vertical GNSS velocities are in a stable northern neighbors reference frame (Liang et al., 2013), both transformed from velocities in the IRTF2008 reference frame (Text S4). For both ascending and descending tracks, the heading angle, $\phi$, is defined as positive clockwise from north, and the incidence angle, $\theta$, is defined as the angle between the LOS and the vertical (Figure S7 in Supporting Information S1).

To strike a balance between quality and abundance of GNSS data points, we used only 3D stations (with east, north and vertical components, Figure S5) when there are more than 10 of them in a frame. For three frames with fewer 3D stations, we included 2D stations (with only east and north components, Figure S5), assuming zero vertical velocities ($V_U = 0$) at these stations. Given the near zero mean of the contribution to LOS from $V_U$ estimated at the 3D sites and their almost equivalent-in-magnitude uncertainties, this is a reasonable assumption. To reflect the lack of information on $V_U$ at the 2D stations, we set the uncertainty in the vertical component, $\sigma_U$, to be the largest absolute value of the available $V_U$ in the frame. As such, the magnitude of uncertainties associated with the GNSS LOS are comparable to that of the InSAR LOS, both at the ~1 mm/yr level. This means the combined uncertainties of InSAR-InSAR offsets and that of InSAR-GNSS offsets are also of similar magnitudes. Therefore, we only needed to prevent the joint inversion from being dominated by the much greater number of InSAR-InSAR offsets relative to the InSAR-GNSS offsets ($g \gg h$). To do so, we weight the bottom three lines (i.e., lines involving GNSS) in Equation 2 by $\sqrt{g/h}$, so that the sum of the mean-squared misfits from these two data sets is minimized.

## 3.2. Effectiveness of Joint Inversion

After mosaicing the InSAR LOS frames using the method described above, the stitched InSAR LOS tracks no longer show along-track frame boundaries (Figure 3). Instead, the track boundaries become more obvious, reflecting the across-track change in incidence angle. From north to south, the ascending LOS velocities become increasingly negative and the descending LOS velocities become increasingly positive, as expected from the increasing east-component velocities (Figure 1a). The remaining LOS changes are well-aligned with mapped faults, especially along the single-strand portions of the Haiyuan and East Kunlun Faults.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Ascending and (b) descending line-of-sight (LOS) velocity tracks in the reference frame provided by the GNSS velocities (Liang et al., 2013; M. Wang & Shen, 2020). GNSS LOS are shown in colored dots.</figcaption></figure>

Statistically, the resultant along-track frame overlap residuals have standard deviations that range between 0.3 and 2.1 mm/yr and average to 0.87 mm/yr (Figure S8 in Supporting Information S1). Overlaps with larger residuals are associated with frames impacted by desert in the northeast, and vegetation and loess in the southeast. The histograms of these frame overlap residuals also follow Gaussian distributions centered near zero (average mode of -0.02 mm/yr), characteristic of random noise. The InSAR-GNSS LOS residuals, with on average 45 GNSS points per track, have an average mode of 0.04 mm/yr and an average standard deviation of 1.27 mm/yr (Figure S9 in Supporting Information S1). If equally split between the InSAR and GNSS, this InSAR-GNSS LOS residual suggests less than 1 mm/yr error in each data set.

After the joint inversion, the across-track LOS differences, subtracting each eastern track from the neighboring western track, show negative modes averaging -1.4 mm/yr (Figure S8 in Supporting Information S1). The systematic offsets are consistent with the change in incidence angles across track. The offsets in the track overlaps also have higher standard deviations (averaging 1.6 mm/yr) than the frame overlaps because of the completely independent interferogram networks and atmospheric delays between tracks. The histograms of the track overlap offsets are less Gaussian in shape because the differences in incidence angles are scaled by $V_E$ and $V_U$, which are variable along track.

## 4. Velocity Decomposition

The ascending and descending LOS velocity maps provide measurements of velocity in two directions. In order to decompose the two LOS velocities into standard 3-component velocities (east, north and vertical), a third constraint must be provided. Hussain et al. (2018) and Weiss et al. (2020) used the interpolated north velocities from GNSS data and the two LOS velocities to simultaneously solve for the east ($V_E$) and vertical ($V_U$) velocity components. Tymofeyeva and Fialko (2018) and Xu et al. (2021) decomposed the LOS velocities into a vertical component and a horizontal component in the local direction determined by interpolated east and north GNSS velocity fields. Both methods require a good distribution of GNSS velocities, which is not available everywhere on the continents.

L. Shen (2020) and L. Shen et al. (2021) developed an alternative method to obtain $V_E$ straight away from InSAR LOS by decomposing the ascending and descending LOS into $V_E$ and $V_{UN}$, where $V_{UN}$ is the combined LOS contributions from the $V_N$ and $V_U$ components projected onto the north-up plane. This method is a simple way of getting high resolution decomposed fields without using an interpolated GNSS $V_N$. However, it also introduces a different source of error as the $V_{UN}$ vector direction changes slightly with the line-of-sight, the decomposed $V_{UN}$
becoming an approximation of the true values on different tracks. The error from the approximation of $V_{UN}$ can also spill into the $V_E$ obtained. We implement this method and show that the error introduced by the $V_{UN}$ approximation is negligible in comparison to the noise level of the data and build upon this method to extract the $V_U$ from $V_{UN}$ using interpolated $V_N$.

## 4.1. Stage 1: Decomposition Into $V_E$ and $V_{UN}$

Following L. Shen (2020) and L. Shen et al. (2021), we first decomposed the two LOS velocities into $V_E$ and $V_{UN}$ via:

$$V(LOS) = -V_E \cos(\phi) \sin(\theta) + V_{UN} \sqrt{1 - \sin^2(\theta) \cos^2(\phi)} \quad (5)$$

based on the mathematical derivation provided in Figure S7. This partitioning is performed independently for each pixel covered by at least one ascending track and one descending track. The resultant $V_E$ and $V_{UN}$ maps (Figure S10a and S10c in Supporting Information S1) have the shapes of the overlapping area of the ascending and descending data coverage (Figure 1a). The general form of the inversion is as follows:

$$\begin{bmatrix} \frac{1}{\sigma(a_0)} LOS_{a_0} \\ \frac{1}{\sigma(a_1)} LOS_{a_1} \\ \frac{1}{\sigma(d_0)} LOS_{d_0} \\ \frac{1}{\sigma(d_1)} LOS_{d_1} \end{bmatrix} = \begin{bmatrix} \frac{1}{\sigma(a_0)} (-\cos(\phi_{a_0}) \sin(\theta_{a_0}) & \sqrt{1 - \sin^2(\theta_{a_0}) \cos^2(\phi_{a_0})} \\ \frac{1}{\sigma(a_1)} (-\cos(\phi_{a_1}) \sin(\theta_{a_1}) & \sqrt{1 - \sin^2(\theta_{a_1}) \cos^2(\phi_{a_1})} \\ \frac{1}{\sigma(d_0)} (-\cos(\phi_{d_0}) \sin(\theta_{d_0}) & \sqrt{1 - \sin^2(\theta_{d_0}) \cos^2(\phi_{d_0})} \\ \frac{1}{\sigma(d_1)} (-\cos(\phi_{d_1}) \sin(\theta_{d_1}) & \sqrt{1 - \sin^2(\theta_{d_1}) \cos^2(\phi_{d_1})} \end{bmatrix} \begin{bmatrix} V_E \\ V_{UN} \end{bmatrix} \quad (6)$$

where $a_0$, $a_1$, $b_0$ and $b_1$ are notations for the up-to-two possible overlapping ascending tracks and up-to-two possible overlapping descending tracks covering a pixel. $\sigma$ is $\sigma(LOS)$ scaled on a spherical model fit to the $\sigma(LOS)$ profile away from the reference center (Text S5, Figure S11 in Supporting Information S1), such that the decomposition will not be overly influenced by data that happen to be in or near the reference window. The number of rows in Equation 6 per pixel depends on the number of tracks covering the pixel with valid data. By performing the velocity decomposition pixel by pixel instead of across the whole area with a covariance matrix estimated from a semivariogram (i.e., Sudhaus & Jónsson, 2009), we make the computation practical for the large data set.

We propagate the uncertainties $\sigma(LOS)$ to $\sigma(V_E)$ and $\sigma(V_{UN})$ (Figure S10b,S10d in Supporting Information S1) through a data covariance matrix, $cov(LOS)$, with only squares of the scaled $\sigma(LOS)$ (Text S5, Figure S11 in Supporting Information S1) in the diagonal terms, based on the assumption that all LOS measurements (from different tracks) are independent. $\sigma(V_E)$ and $\sigma(V_{UN})$ are obtained from the square-root of the diagonal terms of the covariance matrix of $V_E$ and $V_{UN}$, $cov(V_E, V_{UN}) = [A^T cov(LOS)^{-1} A]^{-1}$, with A the transformation matrix from Equation 6.

This method is suitable for areas away from the poles where the heading direction of the track is near north-south, hence the variation of the angle, $\gamma$, between the $V_{UN}$ vector and the vertical axis is small (Figure S7 in Supporting Information S1). In our study area, $\gamma$ varies between ~6–9° which translates into errors of $V_{UN}$ as the vector sum of ~2% of $V_N$ and ~0.5% of $V_U$, as $V_{UN} = -V_N \sin(\gamma) + V_U \cos(\gamma)$. Therefore, the majority of the error actually comes from the north-south motion. For LOS data of 1 mm/yr uncertainty, this method is robust up to $V_N = 50$ mm/yr.

We further compare our $V_E$ (Figure 4a) with that obtained through the conventional way where the interpolated GNSS $V_N$ is used to remove the contribution of $V_N$ from the two LOS velocities, which are then used to estimate $V_E$ and $V_U$ (Hussain et al., 2018; Weiss et al., 2020). The differences in the two $V_E$ maps are small (0.01 ± 0.03 mm/yr) (Text S6, Figure S12 in Supporting Information S1), but the advantages of the $V_E - V_{UN}$ decomposition are that the results are not biased from, and the resolution not reduced because of, the GNSS interpolation.

As the uncertainty level continues to decrease, with longer time series enabled by Sentinel-1's 20 years mission, the errors introduced by the $V_E$ and $V_{UN}$ decomposition will become more significant. One way to mitigate this effect is to choose a reference frame in which $V_N$ is small. This can be done by transforming the GNSS velocities
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) East velocity, $V_E$, map from direct decomposition of ascending and descending line-of-sight velocity tracks. (b) Vertical velocity, $V_U$, map decomposed from Interferometric Synthetic Aperture Radar (InSAR) $V_{UN}$ (Figure S10b) using interpolated GNSS $V_N$ (Figure S13a). See uncertainty maps in Figure S10 and Figure S13. Elevation is from Shuttle Radar Topography Mission 3-arcsecond elevation model (Farr et al., 2007). Faults: HYF = Haiyuan, QLS = Qilianshan, EKLF = East Kunlun, WQF = West Qinling, GLF = Gulang, ZLHF = Zhuanglanghe, MXSF = Maxianshan, JS = Jishishan, LJ = Lajishan, RYSF = Riyueshan, DBSF = Dabanshan, QHNS = Qinghai Nan Shan, ELSF = Elashan, BLJF = Bailongjiang, TZF = Tazang, AWCF = Awancang, HNF = Hanan, YMSF = Yumushan. Cities: JQ = Jiuquan, ZY = Zhangye, JC = Jinchang, WW = Wuwei, GD = Guide, MQ = Maqu, QHL = Qinghai Lake (cf) Regressions and histograms of differences between InSAR and GNSS $V_E$ and $V_U$ velocities (c and e) Crosses = InSAR and GNSS velocities with uncertainties. Orange line = weighted best-fit linear model. Blue dotted line = one-to-one line for reference (d and f) Differences are calculated as GNSS-InSAR. m = mean, s = standard deviation, n = number of GNSS component velocities in InSAR coverage.</figcaption></figure>

to a reference at the center of the study region (L. Shen, 2020; L. Shen et al., 2021) before using the GNSS to mosaic the InSAR LOS (Section 3).

## 4.2. Stage 2: Decomposition of $V_{UN}$ Into $V_U$ and $V_N$

By this stage, we have isolated $V_E$ (Figure 4a) and are left with the subvertical component $V_{UN}$ (Figure S10c in Supporting Information S1) which is related to the two as yet undetermined velocity components $V_U$ and $V_N$. It is now necessary to introduce external constraints on the velocities, that is, GNSS, to separate these two components. As InSAR is more sensitive to $V_U$ than to $V_N$ and GNSS have higher quality measurements and more data points in $V_N$ than in $V_U$, we interpolated $V_N$ from GNSS velocities to obtain the $V_U$ field.

We interpolate the north component of the GNSS velocities from (M. Wang & Shen, 2020) using the universal kriging algorithm of the Pykrige (v1.6.0) Python module (Murphy et al., 2021). We first cleaned the GNSS data by removing data points with north uncertainties greater than 0.7 mm/yr, which reduced the majority of cases where nearby points had contradictory measurements. The universal kriging algorithm allows a polynomial surface to be removed before and added back after interpolation. We chose to remove a best-fit third order polynomial surface that best describes the distribution of GNSS $V_N$ values. To determine the distance-dependent weighting of data points for interpolation, we fitted a spherical model to the semivariogram calculated from the detrended $V_N$ values:
$$ \begin{cases} n + p \times \left( \frac{3d}{2r} - \frac{d^3}{2r^3} \right) & d \le r \\ n + p & d > r \end{cases} \qquad (7) $$

where n is the nugget, s is the sill, $p = s - n$ is the partial sill, r is the range, and d is the distance between any pair of GNSS points. A spherical model was chosen instead of a Gaussian or an exponential model to strike a balance between interpolation smoothness and feature retention. This model also determines how uncertainties increase with distance from the data points (Figure S13a and S13b in Supporting Information S1).

For each track that covers a pixel, a velocity $V_U$ and its uncertainty $\sigma(V_U)$ can be calculated as:

$$ V_U = V_{UN} \frac{\sqrt{1 - \sin^2(\theta) \cos^2(\phi)}}{\cos(\theta)} - V_N \sin(\phi) \tan(\theta) \qquad (8) $$

$$ \sigma(V_U) = \sqrt{\sigma^2(V_{UN}) \frac{1 - \sin^2(\theta) \cos^2(\phi)}{\cos^2(\theta)} + \sigma^2(V_N) \frac{\sin^2(\phi) \sin^2(\theta)}{\cos^2(\theta)}} \qquad (9) $$

hence, pixels in the track overlaps can have up to four realizations of $V_U$ and $\sigma(V_U)$. For such pixels, a weighted mean of the $V_U$, $\bar{V}_U = \frac{\sum V_U i / \sigma(V_U)_i}{\sum 1 / \sigma(V_U)_i^2}$, and the uncertainty of the weighted mean, $\sigma(\bar{V}_U) = 1 / \sqrt{\sum 1 / \sigma(V_U)_i^2}$, are evaluated for the final result (Figure 4b and Figure S13c and S13d in Supporting Information S1).

### 4.3. Comparing InSAR and GNSS Velocities and Uncertainties

As only the LOS components from 3D GNSS stations were used in the frame stitching, except for three frames on the edges with fewer than 10 3D stations (Figure S5), and only the $V_N$ GNSS velocities were used in the second step of velocity decomposition, as an additional quality check, it is possible to compare directly the resultant $V_E$, $V_U$ and their associated uncertainties derived from InSAR with those of the GNSS. The one-to-one fit between InSAR and all 201 stations of GNSS $V_E$ has a gradient of 0.99 and a $R^2$ of 0.84 (Figure 4c) (or a gradient of 0.93 and a $R^2$ of 0.76 if only comparing with 62 2D GNSS stations that were not used for tying InSAR LOS). This agreement is remarkable given the high precision of GNSS velocities in the east component (Figure S10a and S10b, S10e and S10f). In contrast, the $V_U$ components from InSAR and GNSS only exhibit a weak positive correlation with a gradient of 0.33 and $R^2$ of 0.11 (Figure 4e). This fit is not surprising given that the uncertainties associated with the GNSS $V_U$ are almost as large as their absolute values. The better agreement between InSAR and GNSS on the $V_E$ than that on the $V_U$ components can be visualized along two 40 km-wide 600 km-long profiles (Figures S10e, S10f, S13e and S13f). Larger misfits between InSAR and GNSS $V_U$ are found where the GNSS $V_U$ have larger uncertainties, as intended with our weighting strategy in mosaicing the InSAR LOS frames (Section 3.1).

It is worth noting that the InSAR $\sigma V_U$ is only one third of GNSS $\sigma V_U$ (Figure S13d), thanks to the steep line-of-sight of the satellites. The InSAR $V_U$ map is also characterized by a high degree of spatial variability, which would not have been captured by the sparse GNSS points alone, highlighting the advantages of using InSAR for observing vertical motions. The average differences between InSAR and GNSS $V_E$ and $V_U$ are 0.02 mm/yr and -0.05 mm/yr, respectively (Figures 4d and 4f), confirming that our referencing procedure does not add any biases. The standard deviations of the $V_E$ and $V_U$ differences are 1.14 mm/yr and 1.38 mm/yr, respectively, similar to the 1.00 ± 0.26 mm/yr and 1.67 ± 0.23 mm/yr expected from the combined uncertainties from InSAR and GNSS, suggesting the uncertainties estimated for InSAR $V_E$ and $V_U$ are reasonable.

## 5. Strain Rate Calculation

Most existing methods for strain rate calculation were designed to derive continuous strain rate fields from sparse and irregularly spaced leveling or GNSS velocity points (e.g., Z. K. Shen et al., 1996). This problem is typically treated either as a velocity interpolation problem with a user choice of smoothing, coupling or tension
factor (e.g., Sandwell & Wessel, 2016; Smith & Wessel, 1990), or an inversion problem where strain rates or stresses are directly inverted from velocity point data (e.g., Haines et al., 2015; Haines & Holt, 1993; Haines & Wallace, 2020; Kreemer et al., 2014). The former tends to suffer from poorly constrained uncertainties whereas the latter can be limited by computational power. Kinematic block modeling (e.g., Loveless & Meade, 2010; McCaffrey et al., 2013; Meade & Hager, 2005) is yet another approach that is often used to focus strain on predefined faults although the number of blocks tends to grow with the number of velocity measurements incorporated into the model. Pagani et al. (2021) developed a Bayesian approach to extract strain rate from an ensemble of interpolated velocities derived from variable triangular meshes. However, these approaches all rely on some degree of low-pass spatial filtering, hence do not incorporate the full information content of InSAR.

Hitherto, strain rate calculations from InSAR velocities have relied on downsampling the InSAR data because of the computational limit of the N data by N data matrix inversion and the need to reduce short-wavelength noise in the InSAR velocities (Song et al., 2019; H. Wang et al., 2019; Weiss et al., 2020; Xu et al., 2021). Xu et al. (2021) subsampled the decomposed horizontal InSAR velocities to 2.5 km resolution before interpolating with the *gpsgridder* algorithm of Sandwell and Wessel (2016). For studies using the VELMAP method (H. Wang & Wright, 2012), strain rate is obtained from velocities inverted at the vertices of a triangular mesh using spherical approximation equations (Savage et al., 2001). As the inversion is underdetermined, it needs to be regularized through a smoothing criteria. Although the strength of smoothing is typically chosen by examining the trade-off between misfit and solution roughness, the exact choice of smoothing strength remains somewhat subjective. The data downsampling, smoothing, and low-resolution mesh act to dampen the true strain rate magnitudes and mask large velocities discontinuities that might be associated with creeping faults. Here, we present a new method for strain calculation that preserves the resolution of the InSAR velocities and the sharpness of the strain patterns as much as possible and guides noise suppression using synthetic tests performed on 1D velocity profiles.

## 5.1. Combining Gradients of Filtered $V_E$ and Interpolated $V_N$

The gradients of the unfiltered velocities are very noisy due to the ~1 mm/yr uncertainty level of the pixels and the close spacing between velocity observations, thus various filtering strategies were tested. To investigate how filtering affects the signal-to-noise ratio and the strain localisation on faults, we performed synthetic tests on a velocity profile generated using a 1D elastic dislocation model (i.e., a screw-dislocation; Savage & Burford, 1973) at the same resolution, and with a similar noise level, as the InSAR $V_E$ estimates (Figures S14 and S15). Results show that a sliding median filter is better at preserving the shapes of the velocity and gradient profiles than a mean filter (Figure S14a). In addition, applying the filter to the velocity profile before taking the gradients produces more stable results than applying the filter to the gradients of the unfiltered velocity profile (Figure S14b). Subsampling the ~100 m resolution velocity profile to ~1 km resolution does not alter profile shape or affect inverted parameters but effectively reduces computational time.

We also find that the filter window required to suppress noise is a function of the velocity step, fault locking depth and noise level of the data. Noise in the velocity profile can result in velocity gradients of oscillating signs, which when squared, can result in overestimation of the maximum shear and second invariant of strain rates (Equation 10). In order to preserve the total absolute gradients along an arctangent velocity profile, a larger filter window is needed not only for a profile with a higher noise level, but also for one with a smaller velocity step or a shallower fault locking, as they both produce smaller velocity gradients along longer stretches of the synthetic profile, which can be dwarfed by noise (Figures S15 and S16). We experimented with different window sizes and found that a 60 km window is the most suitable for smoothing velocity profiles with noise between 0.5 and 1 mm/yr, velocity steps between 3 and 15 mm/yr, typical of slip rates across the Haiyuan and East Kunlun Faults (Kirby et al., 2007; C. Li et al., 2009; J. Li et al., 2016; Yao et al., 2019; Shao et al., 2020), and locking depths between 3 and 20, typical of locking depths for continental crust (Wright et al., 2013) (Figures S15 and S16 in Supporting Information S1).

Therefore, we applied a sliding median filter with a 60 × 60 km window size to the InSAR $V_E$ (Figure S17) before calculating the gradients $\partial V_E/\partial x$ and $\partial V_E/\partial y$. Because the interpolated GNSS $V_N$ is already smooth due to the sparse distribution of control points (Figure S13a,b), we directly evaluate $\partial V_N/\partial x$ and $\partial V_N/\partial y$ without applying additional filtering. The horizontal strain-rate tensor, $\dot{\epsilon}_h$, is the symmetrical component of the velocity gradient tensor composed of the four velocity gradients, so that the horizontal dilatation rate ($\dot{\epsilon}_{dil}$, Figure 5e), the maximum
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) $\partial V_N/\partial y$ and (b) $\partial V_N/\partial x$ velocity gradients calculated from the interpolated GNSS $V_N$ map shown in Figure S13a. (c) $\partial V_E/\partial y$ and (d) $\partial V_E/\partial x$ calculated from the $V_E$ map median-filtered with a 60 × 60 km window size. Maps of (e) dilatation rate, (f) maximum shear and (g) second invariant of the horizontal strain-rate tensor, calculated from the velocity gradients (ad) using Equation 10. The color scales are chosen to highlight strain concentration along faults.</figcaption></figure>

shear rate ($\dot{\epsilon}_{shear}$, Figure 5f), and the second invariant of the horizontal strain-rate tensor ($\dot{\epsilon}_{IIh}$, Figure 5g), are
defined as follows (Turcotte & Schubert, 2014; Sandwell & Wessel, 2016; M. Wang & Shen, 2020):

$$ \dot{\epsilon}_h = \begin{bmatrix} \dot{\epsilon}_{xx} & \dot{\epsilon}_{xy} \\ \dot{\epsilon}_{yx} & \dot{\epsilon}_{yy} \end{bmatrix} = \begin{bmatrix} \frac{\partial V_E}{\partial x} & \frac{1}{2} \left( \frac{\partial V_E}{\partial y} + \frac{\partial V_N}{\partial x} \right) \\ \frac{1}{2} \left( \frac{\partial V_N}{\partial x} + \frac{\partial V_E}{\partial y} \right) & \frac{\partial V_N}{\partial y} \end{bmatrix} \quad (10) $$

$$ \dot{\epsilon}_{shear} = \sqrt{\dot{\epsilon}_{xy}^2 + (\dot{\epsilon}_{xx} - \dot{\epsilon}_{yy})^2/4} $$

$$ \dot{\epsilon}_{IIh} = \sqrt{\dot{\epsilon}_{xx}^2 + 2\dot{\epsilon}_{xy}^2 + \dot{\epsilon}_{yy}^2} $$

## 5.2. Uncertainties in the Strain Rates

### 5.2.1. Quantitative Estimates of Strain Rate Uncertainties

The uncertainties in the strain rates should come from that of the component velocity gradients, which in turn depend on the uncertainties of the median filtered InSAR $V_E$ and interpolated GNSS $V_N$ fields. Analytically deriving the uncertainties of $V_E$ gradients taking into account the spatial covariance between pixels is challenging because of the median filtering step. Therefore, we empirically estimated $V_E$ gradient uncertainties at an area overlapped by two ascending and two descending tracks. We first created four versions of $V_E$, via $V_E-V_{UN}$ decomposition (Section 4.1), from four possible combinations of one of the ascending tracks and one of the
descending tracks. Then, we median-filtered them and compared their horizontal gradients (Figure S18). The four versions of $V_E$ gradients give standard deviations of 7.8 nst/yr for $\partial V_E/\partial x$ and 4.3 nst/yr for $\partial V_E/\partial y$, which represent 17% and 10% uncertainty levels relative to the gradient signals (46.4 nst/yr for $\partial V_E/\partial x$ and 45.9 nst/yr for $\partial V_E/\partial y$) (Figure S18). These uncertainties are upper bounds for two reasons. First, the error introduced by the approximation of the $V_{UN}$ vector is the largest on the edges, as in the track overlaps, as compared to in the central parts of a track. Second, our $V_E$ (Figure 4a) in track overlaps is evaluated from all available (three or four) tracks instead of just two tracks used in this experiment; the addition of neighboring track(s) would reduce this error with both near-range and far-range LOS velocities involved in the inversion (Section 4.1, Figure S12 in Supporting Information S1).

To estimate the uncertainties of $V_N$ gradients, we used the Monte Carlo method to create 100 interpolated GNSS $V_N$ fields, each time varying each GNSS $V_N$ velocity by an amount randomly sampled from a Gaussian distribution with associated $\sigma(V_N)$, assuming the $\sigma(V_N)$ are independent. From the 100 derived $\partial V_N/\partial x$ and $\partial V_N/\partial y$ fields, we find 7.2 nst/yr and 5.7 nst/yr of uncertainties relative to 10.6 nst/yr and 10.8 nst/yr of the gradient signals, respectively (Figure S19). These statistics reflect the high noise-to-signal ratios (68% for $\partial V_N/\partial x$ and 53% for $\partial V_N/\partial y$) of the gradients of a spatially smooth $V_N$ field. Taking 7.8 nst/yr for $\sigma(\partial V_E/\partial x)$, 4.3 nst/yr for $\sigma(\partial V_E/\partial y)$, 7.2 nst/yr for $\sigma(\partial V_N/\partial x)$ and 5.7 nst/yr for $\sigma(\partial V_N/\partial y)$, we have uncertainties of 6.3 nst/yr on the maximum shear rate, 9.7 nst/yr on the dilatation rate and 7.0 nst/yr on the second invariant.

### 5.2.2. Overestimated Strain in Areas of Large $\sigma(V_E)$

As the velocity gradients are squared in the evaluation of maximum shear and second invariant, noise remaining in the filtered $V_E$ will lead to oscillating gradients of opposite signs (Figure S15) and an overestimation of the strain rates. This is most likely the case around the West Qinling Faults and near the eastern termination of the East Kunlun Fault (Figure 5), where $\sigma(V_E)$ is up to 2 mm/yr (Figure S10b). We estimate the degree of this overestimation by filtering 1000 simulations of arctangent velocity profiles and analyzing the absolute $V_E$ gradients summed over the filtered profile relative to that of the model. For a profile with 2 mm/yr Gaussian noise, 20 km locking depth and 5 mm/yr of velocity step, using a 60 km filter window can lead to overestimates in the absolute gradients by $23 \pm 8\%$. This percentage overestimation will increase with smaller velocity steps and shallower locking depths as the signal-to-noise ratio decreases.

### 5.2.3. Underestimated Strain From Low-Resolution GNSS $V_N$ Gradients

The interpolated GNSS $V_N$ field is smoother than the $V_E$ field due to the lower density of the GNSS measurements. As a result, the GNSS velocity gradients are also less sharp and have lower amplitudes than the $V_E$ gradients (Figures 5a–5d). This difference has an impact on the spatial patterns in the strain-rate maps, with contributions from the $V_N$ gradients being comparatively weaker. Therefore, we should be mindful of potentially underestimated shear strain rates over N-S oriented strike-slip faults (e.g., Riyueshan Fault) and contractional strain rates over E-W oriented reverse faults (e.g., Qilianshan thrusts).

### 5.2.4. Location Uncertainty From Uneven Pixel Densities

When the pixel density of the InSAR $V_E$ is not uniform, the median filter may result in lateral shifts of the strain rate features toward the side with poorer pixel coverage by distances up to half the filter window size (30 km in our case). Pixel decorrelation between areas of different velocities may also alter the shapes of the strain rate features. However, the total strain rate across profiles and across the mapping area should be conserved, and no strain rate peak should appear without a change in velocity.

## 6. Main Features in the Velocity and Strain-Rate Fields

Our analysis results in two regional velocity maps (Figures 4a and 4b) covering an area of 440,000 km² with ~100 m resolution, as well as three strain rate maps at ~1 km resolution with better spatial coverage because of the filtering (Figures 5e–5g). These maps are consistent with previously published GNSS-based strain rate products (M. Wang & Shen, 2020) in the long wavelength (Figure S20), but provide a new, extensive view of the crustal motion of the NE Tibetan Plateau at unprecedented resolution.
## 6.1. Strain Localisation Along Haiyuan and East Kunlun Faults

We observe concentrated strain along the major Haiyuan and East Kunlun Faults (Figure 5g). We also see more distributed strain across the Qilianshan and West Qinling Fault, as well as some finite strain away from the main faults. If we mask out pixels with second invariant strain rate above 40 nst/yr, as highlighted in Figure 5g, the remaining pixels together sum up to half of the total strain rate in the study region and average to 21 nst/yr. The major faults accommodate on average 76 nst/yr of strain rate across 20% of the area.

The spatial extent of the strike-slip motion on the Haiyuan Fault is most clearly defined in the maximum shear rate map (Figure 5f), with high shear extending from Haiyuan County in the east, to Qilian County in the west, where the strike-slip fault terminates and transitions into the Changma and North Qilianshan thrusts. This is consistent with the change of strain patterns to the west of Qilian County where the maximum shear diminishes and transitions to contraction in the form of strongly negative dilatation rates (Figure 5e). In comparison to the Haiyuan Fault, elevated strain on the East Kunlun Fault is less tightly localized on the main fault strand, with shear observed to the south of the fault near Maqing and to the north of the fault near Maqu, where the strain appears distributed across multiple strands.

## 6.2. Contractional Strain at Fault Junctions

We also observe focused strain at fault junctions where the Haiyuan Fault branches into the Qilianshan Thrusts and where it interacts with the Gulang and Lianhuashan Faults, as well as where the East Kunlun Fault meets the Awangcang Fault. These transition areas are characterized by contractional strain rates (i.e., negative dilatation; Figure 5e), suggesting volumetric changes and stress anomalies at the triple junctions, likely induced by relative motion of neighboring crustal material at the branching points. This observation is in agreement with the model of King and Nábělek (1985) where off-fault strain is expected at a process zone around a three-fault junction to permit zero displacement at the fault intersection. We cannot tell the difference between elastic strain accumulation and plastic strain that does not end up in earthquakes. Further modeling could assess the degree to which off fault deformation is required to release the strain observed. Understanding such focused strain at the fault bends and branches is crucial as such geometric complexities might promote the initiation and termination of future earthquakes (Bhat et al., 2007; Quigley et al., 2019; Sathiakumar & Barbot, 2021; Schwartz & Sibson, 1989; Walters et al., 2018; Wesnousky, 2008).

## 6.3. Strain-Rate and $V_E$ Profiles Across Faults

We further compare strain rates and $V_E$ across the Qileshan and Haiyuan, East Kunlun and West Qinling Faults using a suite of fault-perpendicular profiles (Figure 6). The magnitude and shapes of the strain peaks are largely related to the step size and gradients of the $V_E$ profiles. A typical arctangent-shaped profile over the East Kunlun Fault (Figure 6c, profile 2-2') corresponds with a Gaussian-shaped strain-rate profile. The 7.3 mm/yr velocity step is likely a combined contribution from both East Kunlun and Awangcang Faults. Despite having a larger velocity step, this long wavelength $V_E$ profile gives rise to a broad strain profile with a lower peak strain value than that of the profile two over the Haiyuan Fault (Figure 6b, profile 2-2'), where strain is more localized. A sharp velocity step in $V_E$ in profile 1 across the Haiyuan Fault corresponds to a spike in the strain-rate profile, with a secondary step over the Gulang Fault indicated by a smaller peak in the strain-rate profile (Figure 6b, profile 1-1'). Multiple strain peaks are observed across both profiles over the Qileshan (Figure 6a), as expected from strain partitioning across the parallel thrusts. The West Qinling Fault is characterized by nearly linear velocity profiles which correspond to near flat strain rate profiles without any prominent peaks associated with the surface fault traces. This suggests that either this area is characterized by distributed deformation, or any localized strain is hidden in the noise caused by vegetation and loess cover in this area, or we could be seeing a seismic cycle effect with rapid postseismic strain rate after a M 8 events that occurred here in 1654 (CEA-EDPB, 1995) and slow interseismic strain rate observed now (Elliott et al., 2016; Hussain et al., 2018; Zhu et al., 2020).

## 6.4. Creep Detection and Creep Rate Measurement Over the Laohushan Section

In Figure 7a, we plot the second invariant of strain rate map (Figure 5g) using an adjusted color map in order to highlight the 30-km-long Laohushan section where creep has been documented over the past three decades
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison between full-resolution $V_E$ and second invariant of horizontal strain rate profiles across (a) the Qilianshan thrusts, (b) the Haiyuan Fault and the Gulang Fault, (c) the East Kunlun Fault, and (d) the West Qinling Fault. Peak strain and the $V_E$ velocity differences between the mean velocity at the ends of the profiles are labeled in each profile panel. Thick and thin black lines on the profiles show median and standard deviations calculated for each 1 km bin. NQLS = North Qilianshan Fault; YMSF = Yumushan Fault; CMF = Changma Fault; SLNSF = Shulunanshan Fault; LLL = Lenglongling Fault; ZLHF = Zhuanglanghe Fault; AWCF = Awancang Fault; MXSF = Maxianshan Fault; LTDCF = Lintan-Dangchang Fault.</figcaption></figure>

(Cavalié et al., 2008; Jolivet et al., 2013, 2015). The $V_E$ profile in Figure 7c shows that the creep has continued between the years 2014 and 2020 with a $V_E$ step of 2.8 ± 0.3 mm/yr, based on the weighted means of the velocities within 5 km on each side of the fault. The equivalent fault parallel creep rate is 2.9 ± 0.3 mm/yr, assuming no vertical or fault perpendicular motions, consistent with the 2.5 ± 0.4 mm/yr obtained by Y. Li et al. (2021) using GNSS, Evisat and leveling data. This is slower than the 5 ± 1 mm/yr measured by Evisat InSAR data between 2004 and 2009 (Jolivet et al., 2012), and the 6.3 ± 2 mm/yr measured by ERS InSAR data between 1993 and 1998 (Cavalié et al., 2008), suggesting recent decrease in the creep velocity. The $V_E$ and $V_U$ decomposition also helps separate horizontal creep from the subsidence signal in the Jingtai Basin. The Jingtai Basin is located at a releasing jog between the Laohushan segment to the west and the 1920 Haiyuan Earthquake surface rupture to the east (Han et al., 2021; Jolivet et al., 2012). The subsidence could be caused by extension across the basin induced by the 1920 Haiyuan Earthquake, which would also have lowered the normal stress and hence promoted creep. The continuous decrease of the creep rate over the past three decades suggests that this creep is transient, possibly in relation with the rupture termination of 1920 Haiyuan Earthquake.

We also identify a ~66 km long high strain segment further west on the Haiyuan Fault. This could be related to a $M_W$ 5.9 Menyuan Earthquake that occurred on 21 January 2016 (Figures 7a and Figure 8d). However, considering that the fault patch that slipped during the earthquake was only 20 km long (Y. Li, Shan, et al., 2016; H. Wang et al., 2017), yet the high-strain segment is three times longer indicates that the earthquake might have induced postseismic creep. The Menyuan Earthquake was a thrust earthquake that occurred on a fault splay attached to the Haiyuan Fault at a depth of 10 km (Y. Li, Jiang, et al., 2016; H. Wang et al., 2017; Y. Zhang et al., 2020). The motion of the pop-up structure might have released normal stress on the Lenglongling Fault thus making stress conditions conducive to creep on the Haiyuan Fault. This postseismic creep might have, in turn, induced a coulomb stress change that promoted the occurrence of the strike-slip $M_{w,6.6}$ 8 Jan 2022 Menyuan Earthquake on the western end of the creeping section. Similarly, the high-strain signal in the southeast corner of Figure 7a is associated with the 2017 $M_s$ 7.0 Jiuzhaigou Earthquake and associated post-earthquake deformation.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) Second invariant of horizontal strain rate. Circles show $M_W > 5$ earthquakes that occurred between 2012 and 2020 from the Global Centroid Moment Tensor Catalog (Dziewonski et al., 1981), with sizes scaled by magnitude. Event year and magnitude with $M_W > 5.5$ are labeled. (b) Zoomed-in view of the strain rate second invariant for the region indicated with the red box in (a), showing the Laohushan creeping section. Panels (c) and (d) are $V_E$ and $V_U$ maps for the same region. (e) shows the A-A' profile from (b). (f) shows the A-A' profile from (c).</figcaption></figure>

## 6.5. Shear Strain on the Eastern Extension of the Dabanshan Fault

Figure 8 shows a close-up view of the region around the central portion of the Haiyuan Fault that highlights elevated maximum shear strain on a fault branch to the south of and parallel to the Lenglongling Fault. The associated velocity step that is clearly visible in both map and profile views indicate relative motion of ~2 mm/yr on the eastern extension of the Dabanshan Fault (Figure 8c). The surface trace of this fault is from China's Seismic Active Fault Survey Data Center and is labeled as a Holocene fault, possibly because of the 1929 M5.5
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Maps over and profiles across the Lenglongling Section of the Haiyuan Fault. (a) $V_E$ median filtered with a window size of 60 km. (b) full-resolution $V_E$, (c) Maximum shear strain rate showing shear on a parallel branch to the south of Lenglongling Fault, (d) Second invariant of strain rate, plotted with a different color scale, showing signal that possibly represents post-seismic creep. (e) Dilatation rate showing contraction signals at fault junctions and focal mechanisms of a 2015 November 22 Mw 5.2 oblique strike-slip earthquake and the January 2016 21 Menyuan Earthquake. (f) Topography over the fault branches around the Lenglongling segment of the Haiyuan Fault. (g) Profile 1-1' and (h) profile 2-2' with full-resolution $V_E$, filtered $V_E$, maximum shear, dilatation and topography data extracted from (a-c, e and f). County names: GL = Gulang, QL = Qilian, MY = Menyuan. Fault names: LHSN = Lianhuashannan, YMS = Yumushan, SH = Sunan, LLL = Lenglongling, DBS = Dabanshan, RYS = Riyueshan. Basin names: MYB = Menyuan Basin, XNB = Xining Basin. Dashed lines in (f) and (h) mark a potential fault suggested by the shear strain, though its location in (f) is only weakly constrained by the sparse $V_E$ pixels around this feature. Alternative fault geometry from (J. Zhang et al., 2009) is plotted in (f) in red for comparison.</figcaption></figure>

historical earthquake that occurred just to the south (https://activefault-datacenter.cn/) (Xu et al., 2016). How and whether this segment connects with the rest of the fault system is not well constrained due to the sparse InSAR data coverage across the mountains (Figure 8f). However, if the Dabanshan and the Haiyuan Faults are connected and capable of rupturing simultaneously, they pose a previously unrecognized and elevated seismic hazard to the region.

## 6.6. Strain Concentration Near the Eastern Termination of the Haiyuan Fault

Additional intriguing features include three high shear strain rate fingers and an arcuate-shaped concentration of contractional strain on the west flank of the Quwu Shan near the eastern termination of the Haiyuan fault (Figure 9a). This feature starts from the east flank of the Hasishan on the Haiyuan Fault and extends from a southern branch (Figure 9b) of the main fault which is commonly believed to be an active fault that did not rupture during the 1920 Haiyuan Earthquake (Cavalié et al., 2008; C. Li et al., 2009; Ren et al., 2016; Matrau et al., 2019; Yao et al., 2019; Ou et al., 2020). Deng et al. (1986) and Gaudemer et al. (1995) suggested the branch extends southeastwards and connects with an oblique reverse fault on the southwestern flank of the Quwu Shan (mapped in Figure 9b) to form a subordinate fault strand in the Haiyuan Fault Zone. However, we do not observe a velocity gradient or elevated strain across the southwestern flank of Quwu Shan. Instead, a $V_E$ gradient and a strong
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Maps and profiles showing contraction and shear branches that stem from the Haiyuan Fault west of the Quwu Shan. (a) Dilatation, (b) maximum shear strain rate, (c) hillshade of 3-arcsecond digital elevation model from Shuttle Radar Topography Mission (Farr et al., 2007) with diverted streams highlighted in cyan, (d) $V_E$, (e) $V_U$ where “Farm” indicates farmland and (f) profiles showing a tilted alluvial fan surface collocated with contractional strain and $V_E$ and $V_U$ gradients.</figcaption></figure>

contraction signal exist along the near north-south trending valley between the Yellow River and the Quwu Shan (Figures 9a, 9d, and 9f). The $V_U$ map also shows a velocity step across the two sides of the valley, with the west uplifting relative to the east, ignoring the subsidence signals along the valley which correspond to the farmlands near villages and are likely caused by water extraction for irrigation (Figures 9e and 9f, see similar subsidence
patterns discussed in Section 6.7.3). This motion might have created the topographic step visible on the hill shade map along the contractional feature (Figure 9c). This is most noticeable near Gaowanxiang, where the western-most 1 km of an otherwise west-dipping coalesced alluvial fan surface (bajada) is tilted to the east (Figure 9f), diverting the local drainage (Figure 9c). We interpret the signal as a growing fold or a young blind west-dipping thrust fault that, together with the series of thrusts and folds oriented sub-parallel to the Liupanshan, helps absorb the termination of the left-lateral motion on the Haiyuan Fault (Duvall & Clark, 2010). The localisation of the contractional strain here also suggests that the Longzhong basin block, which is bounded by the Haiyuan, Liupanshan, West Qinling, Maxianshan and Zhuanglanghe Faults (Figure 4a) (e.g., X. Li et al., 2021), might not be as rigid as previously thought (Y. Wang et al., 2017; W. Wang et al., 2017).

## 6.7. Vertical Motions From Climatic, Hydrological and Anthropogenic Processes

Unlike the $V_E$ map, which undoubtedly primarily reflects active tectonics, we do not see much correspondence between vertical deformation patterns and the surface fault geometry (Figures 4a and 4b). Instead, we observe a first order correlation between vertical motion and topography. In the high plateau, the $V_U$ maps is dominated by subsidence. In the lowland, the $V_U$ signal is dominated by uplift. The short-wavelength variations in vertical motion compared to $V_E$ suggests the source of the vertical motions are likely shallow and mostly driven by surface or near-surface processes. By comparing $V_U$ with Esri's 2020 land cover map (Figure S21) (Karra et al., 2021), we find strong correlations between $V_U$ and three hydrology-related land cover classes, despite the low percentage of these pixels compared to the overall map. Water (mostly along river courses) and flooded vegetation generally correspond to uplift ($V_U = 1.3 \pm 2.2$ mm/yr), and snow/ice corresponds to subsidence ($V_U = -2.5 \pm 2.9$ mm/yr).

### 6.7.1. Subsidence of Permafrost-Rich Qilianshan

Pixels classified as snow/ice are mostly in the peaks of Qilianshan (Figure S21 in Supporting Information S1), east of where Daout et al. (2020) observed active-layer freeze and thaw cycles and widespread subsidence interpreted as ice-loss. We compare our $V_U$ map with the permafrost zonation index (PZI) map from Gruber (2012), which is an estimate of the degree to which permafrost exists in a region. High PZI means permafrost is expected nearly everywhere, whereas low PZI means permafrost likely exists only in the most favorable conditions (Gruber, 2012). We find that the area of strong subsidence in the Qilianshan, where $V_U < -2$ mm/yr, matches the area where PZI >0.5 (Figures 10a and 10c), suggesting that isotropic thaw subsidence might also span this region of the Tibetan Plateau (Zou et al., 2017; Bibi et al., 2018; Cao et al., 2019; Gao et al., 2021; J. Wang et al., 2021; J. Chen et al., 2022).

### 6.7.2. Uplift From Blocked Drainage

Significant local uplift is detected to the east of Jiuquan, to the north of Jinchang, to the east of Maqu and to the east of Guide cities (Figure 4b). If we compare regions where $V_U > 1.8$ mm/yr to the topography, we find that the rapidly uplifting regions correspond to locations where drainage is blocked by topographic barriers (Figure 10d,e,g,i). The patch near Jiuquan is drained by the Beida and Black Rivers as well as the numerous alluvial fans on the northeast slope of the Qilianshan (Figure 10d). The outlet is blocked by the Jintanan Shan and Heli Shan in the north such that the water is ponding in front of the linear ridge. Similarly, the triangular patch to the north of Jinchang sits in the Chaoshui Basin, where an endorheic river drained from the northwest terminates. This region is sandwiched between the alluvial fans stemming from the Beida Shan in the north and Longshou Shan in the south, collecting water from two almost opposite azimuths (Figure 10e). The same is true for the fast uplift rates in the Gonghe Basin (Figures 10a–10c), where the melt water from the permafrost-rich peaks in the west fails to reach the Yellow River and stalls in a local topographic minimum. Further south, high $V_U$ values are observed in the Zoige Basin to the east of Maqu where the White and Black Rivers join the Yellow River and the Yellow River makes a U-turn in front of the Min and East Kunlun Mountains. This is also where pixels classified as flooded vegetation are clustered (Figure 10f), consistent with the wetland environment of the Zoige Basin (B. Li et al., 2014). As the Yellow River enters the Xunhua Basin (Figures 10h and 10i), the outflow is again blocked by the N-S trending Jishi Shan such that the entire eastern half of the Xunhua Basin shows rapid uplift rates. The locations of rapid uplift relative to drainage and topography suggests a clear hydrological origin of the uplifting signals. Further work is required to determine whether this is a manifestation of poroelastic expansion of the soil or the cumulative effect of phase bias due to systematic changes in soil properties. One way or
<figure><img src="image_12.png" /><figcaption>Figure 10.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) $V_U$, (b) $V_E$, (c) dilatation rate and (d) maximum shear rate maps of the Longyangxia Reservoir. (e) Profiles from (a–d) showing subsidence around and convergence toward the center of the reservoir causing dilatation and shear in the area. Thick and thin lines mark binned mean and standard deviations along profiles. Labels on each profile show lowest mean $V_U$, difference in mean $V_E$ on the edges of the profile, lowest mean dilatation and highest mean maximum shear.</figcaption></figure>

another, this phenomena requires a water mass increase likely linked to increases in surface air temperature and precipitation, melting of glaciers and thawing of ice-rich permafrost (e.g., Bibi et al., 2018).

### 6.7.3. Subsidence of Farmlands in the Hexi Corridor

Patches of subsidence are observed along the Hexi Corridor to the north of the permafrost zone (Figures 10a and 10c). The spatial coverage of the subsidence signals corresponds to farmland near the major cities of Jiuquan, Zhangye, Jinchang and Wuwei (Figures 4b and Figure 10c). We attribute the subsidence to rapid groundwater extraction for irrigation which is common practice in the arid northwest of China, where the loss of groundwater storage has long been a known issue (He et al., 2012; X. Zhang et al., 2018; Wu et al., 2019). The fastest subsiding crop areas are those to the east of Jinchang and to the north of Wuwei (Figure 4), where $V_U < -10$ mm/yr are observed.

### 6.7.4. Subsidence Around the Longyangxia Reservoir

The moon-shaped empty patch to the south of the Qinghai Lake is the Longyangxia Reservoir (Figures 10a and 11), and is at the center of a subsiding region. The subsidence signal is strongest at the thin strip of land in the head of the hook-shaped water body, with a subsidence rate of ~5.2 mm/yr (Figure 11e). The subsidence rate decreases away from the reservoir over a wavelength of >20 km on each side, typical of the length scale of elastic deformation of the upper crust (Doin et al., 2015). The $V_E$ profile shows a step that is greater just adjacent

**Figure 10.** Interpretations of $V_U$ signals in (a–c) the Qilianshan, (d and e) the Hexi Corridor, (f and g) the Zoige Basin, and (h, i) the Xunhua Basin. See location of (d) in panel (b) and locations of (e, f and h) in Figure 4b. (a) $V_U$, (b) 3-arcsecond digital elevation model (DEM) from Shuttle Radar Topography Mission, (c) Permafrost Zonation Index (PZI) (Gruber, 2012) map and farmland (crops) in ESRI's global land cover map (Karra et al., 2021). Black lines in (a–c) outline the PZI = 0.5 contour in the Qilianshan and the boundaries of the anomalous uplift and subsidence signals in the Hexi Corridor (d, e, g and i) show pixels with $V_U > 1.8$ mm/yr on top of hill shade with drainage data from Generic Mapping Tools (Wessel et al., 2013) and HydroRIVERS data set (available at https://www.hydrosheds.org) (Lehner & Grill, 2013). (f) shows permanent rivers and flooded vegetation pixels in purple (Karra et al., 2021) on top of topography. (h) shows faults and permanent rivers on top of topography.
to the reservoir and smaller further away, suggesting material is being pulled toward the center of subsidence. Therefore, we interpret subsidence as a result of the gravitational loading from the water column, in response to a ~20 m overall increase in water level detected between 2017 and 2020 (Zhao et al., 2022). This water load might have been established prior to 2017 and caused convergence of the surrounding crust and hence the contraction observed (Figure 11c). The source of the shear strain is unclear (Figures 11c and 11d). One possible explanation is that the interpolated GNSS $V_N$ field is not of high-enough resolution to capture this short-wavelength variation, hence the $V_N$ gradients cannot balance the $V_E$ gradients in the evaluation of the maximum shear rate.

### 6.7.5. Large-Scale Regional Uplift

Focusing only on regions where the absolute value of $V_U$ is less than 2 mm/yr (i.e., ignoring non-tectonic deformation) highlights large-scale, regional uplift northward of the West Qinling Fault, across the Haiyuan Fault, along most of the Hexi Corridor, as well as in the lens-shaped valleys between the Qilian ranges and the low-relief plains to the west of Qinghai Lake (Figure 4b). Whilst the bedrock vertical motions in the high plateaus are still mixed and uncertain, we can see a broad-scale uplifting signal across the low-lying northeast half of the study region (compare Figures 1a and 4b). The uplift is relative to the average $V_U$ velocities of three reference GNSS stations in Ordos Plateau, Gobi-Alashan and Mongolia chosen by Liang et al. (2013) in the stable northern neighbors reference frame. This observation could be related to the northward growth of topography as inferred by the new-found activity on the thrusts in Heli Shan to the north of the Hexi Corridor measured by W. J. Zheng et al. (2013), or could be a result of increased erosion, hence isostatic uplift, due to the increased drainage from climate change.

## 7. Conclusions

In this study, we have derived ~100 m resolution $V_E$ and $V_U$ maps at ~1 mm/yr uncertainty levels and three ~1 km resolution strain rate maps covering 440,000 km² of the NE Tibetan Plateau. We developed new methods for mosaicing frame-based LOS velocity maps into regional velocity maps, deriving Cartesian (i.e., north-south, east-west, and vertical) velocities, and deriving high-resolution strain rate maps from the velocity fields that preserve the rich and detailed information contained in the dense geodetic data sets. The resulting $V_E$ and $V_U$ maps highlight for the first time at a large scale tightly focused left-lateral, strike-slip deformation across the major Haiyuan and East Kunlun fault system, and an outward growth of topography to the north-east of the Tibetan plateau. The second-invariant of the horizontal strain-rate tensor shows that half of the total horizontal strain rate is accumulating on those major structures. Shear strain rate helps delineate the extent of the major strike-slip faults, with a mix of both localized and distributed deformation as a result of fault geometry and slip partitioning. Dilational and compressional strains are observed between Qilianshan thrusts and on fault junctions, such as where the Haiyuan Fault branches into neighboring strike-slip and thrust faults and where the Awancang strike-slip fault merges into the East Kunlun Fault, highlighting the local volumetric changes within zones of complex fault geometry. The high-resolution strain rate maps also aid in identifying fault segments that are creeping at the surface and previously unknown active faults. The Laohushan creeping section of the Haiyuan Fault is mapped to be 30 km long across both the Laohushan with an average creep rate of ~3 mm/yr, which is >20% lower than the creep rate reported in previous studies. High shear-strain rates are also observed along the Lenglongling section of the Haiyuan Fault, potentially induced by postseismic creep following the 2016 Menyuan Earthquake. The eastward motion of the eastern half of the Menyuan Basin also appears to be decoupled from that in the Xining Basin, causing shear strain to accumulate on the eastern extension of Dabanshan Fault. Strong north-south-trending contraction is identified on the western flank of the Quwushan, in agreement with a back-tilted alluvial fan surface. The $V_U$ map reveals rapid subsidence across the permafrost-rich Qilianshan due to large-scale ice loss, at the oases along the Hexi Corridor due to the extraction of groundwater for irrigation, and around the Longyangxia Reservoir due to gravitational loading from the increasing water level. Small-scale uplifts are also observed at the Jiuquan, Chaoshui, Zoige and Xunhua Basins where river drainage hits topographic barriers. With the ongoing expansion of global InSAR data processing and coverage, high-resolution velocity and strain rate maps like those shown here will play an increasingly important role in measuring active fault deformation, detecting unknown and creeping faults, advancing our understanding of the regional and global seismic hazard caused by active faults and characterization of climatic, hydrologic and anthropogenic-related processes alike.
## Data Availability Statement

The original interferograms are available on COMET-LiCS portal (https://comet.nerc.ac.uk/comet-lics-portal/). The line-of-sight and Cartesian velocities, their associated uncertainties, and the strain rate fields that we derived are available at CEDA archive via http://dx.doi.org/10.5285/7fbb95c288af44ab8b40e74fef0e7cbc with open government access. Version 1.0 of the scripts used for mosaicing LOS velocities along track, correcting for the reference effects on LOS uncertainties, interpolating GNSS $V_N$, is preserved at DOI:10.5281/zenodo.6546922, available via the Creative Commons Attribution 4.0 International license.

## Acknowledgments

We would like to thank the European Space Agency for the open data policy on the Sentinel-1 mission. This study was funded by the Natural Environment Research Council (NERC) through the "Looking inside the Contins from Space (LiCS)" large grant (NE/ K011006/1), the Centre for the Obser- vation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET), a partnership between UK Universities and the British Geological Survey, and the "Earthquakes without Frontiers (EwF)" consortium (NE/J01978X/1). Qi Ou has been supported by scholarships and awards from China Scholarship Council-Pacific Alliance Group Oxford Scholarship and Great Britain-China Educational Trust. LiCSAR data contains modified Copernicus Sentinel data (2014–2019) analysed by the COMET-LiCSAR team; we thank Yasser Maghsoudi for his support. LiCSAR processing is carried out at JASMIN, the UK’s collaborative data analysis facility for environmental science (http://jasmin.ac.uk) hosted at the Centre for Environ- mental Data Analysis (CEDA). We would like to thank Paul Tregoning, the editor, Jianbao Sun and another anonymous reviewer for their constructive comments. We thank Min Wang and Shiming Liang for their advice on combining the horizontal and vertical GNSS velocities and thank Min Wang for providing GNSS strain rate data for comparison with our results. We also thank Jérôme Lavé for helpful discussions on the geomorpho- logic observations that corroborated our interpretation of strain features. Figures were plotted using the Generic Mapping Tools Version 6 (Wessel et al., 2019) and Python.

## References

Ansari, H., De Zan, F., & Parizzi, A. (2020). Study of systematic bias in measuring surface deformation with SAR interferometry. *IEEE Transactions on Geoscience and Remote Sensing*, 59(2), 1–12. https://doi.org/10.1109/TGRS.2020.3003421

Bekaert, D. P., Biessel, R., Havazli, E., Hua, H., Karim, M., Lai-Norling, J., et al. (2020). New and upcoming developments of standardized InSAR products by the Advanced Rapid Imaging and Analysis (ARIA) Project for Natural Hazards. *AGU Fall Meeting Abstracts*, 2020, G004–0038.

Bhat, H. S., Olives, M., Dmowska, R., & Rice, J. R. (2007). Role of fault branches in earthquake rupture dynamics. *Journal of Geophysical Research*, 112(11), 1–16. https://doi.org/10.1029/2007JB005027

Bibi, S., Wang, L., Li, X., Zhou, J., Chen, D., & Yao, T. (2018). Climatic and associated cryospheric, biospheric, and hydrological changes on the Tibetan plateau: A review. *International Journal of Climatology*, 38, e1–e17. https://doi.org/10.1002/joc.5411

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

Burgmann, R., Rosen, P. A., & Fielding, E. J. (2000). Synthetic aperture radar interferometry to measure Earth’s surface topography and its deformation. *Annual Review of Earth and Planetary Sciences*, 28(1), 169–209. https://doi.org/10.1146/annurev.earth.28.1.169

Cao, B., Zhang, T., Wu, Q., Sheng, Y., Zhao, L., & Zou, D. (2019). Permafrost zonation index map and statistics over the Qinghai–Tibet Plateau based on field evidence. *Permafrost and Periglacial Processes*, 30(3), 178–194. https://doi.org/10.1002/ppp.2006

Cavalié, O., Lasserre, C., Doin, M. P., Peltzer, G., Sun, J., Xu, X., & Shen, Z. K. (2008). Present-day deformation across the Haiyuan fault (Gansu, China), measured by SAR interferometry. *European Space Agency*, (Special Publication). (655).

Chen, C. W., & Zebker, H. A. (2002). Phase unwrapping for large SAR interferograms: Statistical segmentation and generalized network models. *IEEE Transactions on Geoscience and Remote Sensing*, 40(8), 1709–1719. https://doi.org/10.1109/TGRS.2002.802453

Chen, J., Wu, T., Zou, D., Liu, L., Wu, X., Gong, W., et al. (2022). Magnitudes and patterns of large-scale permafrost ground deformation revealed by Sentinel-1 InSAR on the central Qinghai–Tibet Plateau. *Remote Sensing of Environment*, 268, 112778. https://doi.org/10.1016/j.rse.2021.112778

China Earthquake Administration Earthquake Damage Protection Bureau (CEA-EDPB). (1995). In N. Jiang (Ed.), *China historical strong earthquake catalogue (BC23–AD1911)* (1st ed.). Earthquake Press.

Daout, S., Dini, B., Haeberli, W., Doin, M. P., & Parsons, B. (2020). Ice loss in the Northeastern Tibetan Plateau permafrost as seen by 16 yr of ESA SAR missions. *Earth and Planetary Science Letters*, 545, 116404. https://doi.org/10.1016/j.epsl.2020.116404

Daout, S., Doin, M. P., Peltzer, G., Lasserre, C., Socquet, A., Volat, M., & Sudhaus, H. (2018). Strain partitioning and present-day fault kinematics in NW Tibet from Envisat SAR interferometry. *Journal of Geophysical Research: Solid Earth*, 123(3), 2462–2483. https://doi.org/10.1002/2017JB015020

Daout, S., Doin, M.-P., Peltzer, G., Socquet, A., & Lasserre, C. (2017). Large-scale InSAR monitoring of permafrost freeze-thaw cycles on the Tibetan Plateau. *Geophysical Research Letters*, 44(2), 901–909. https://doi.org/10.1002/2016GL070781

Daout, S., Jolivet, R., Lasserre, C., Doin, M. P., Barbot, S., Tapponnier, P., et al. (2016). Along-strike variations of the partitioning of convergence across the Haiyuan fault system detected by InSAR. *Geophysical Journal International*, 205(1), 536–547. https://doi.org/10.1093/gji/ggw028

Daout, S., Sudhaus, H., Kausch, T., Steinberg, A., & Dini, B. (2019). Interseismic and postseismic shallow creep of the North Qaidam thrust faults detected with a multitemporal InSAR analysis. *Journal of Geophysical Research: Solid Earth*, 124(7), 7259–7279. https://doi.org/10.1029/2019JB017692

Deng, Q., Chen, S., Song, F., Zhu, S., Wang, Y., Zhang, W., et al. (1986). *Variations in the geometry and amount of slip on the Haiyuan (Nanx-ihaushan) Fault zone, China and the surface rupture of the 1920 Haiyuan earthquake*. In *Haiyuan fault zone* (pp. 169–182). https://doi.org/10.1029/GM037p0169

Deng, Q., Sung, F., Zhu, S., Li, M., Wang, T., Zhang, W., et al. (1984). Active faulting and tectonics of the Ningxia-Hui Autonomous Region, China. *Journal of Geophysical Research*, 89(10), 4427–4445. https://doi.org/10.1029/JB089iB06p04427

Doin, M. P., Twardzik, C., Ducret, G., Lasserre, C., Guillaso, S., & Jianbao, S. (2015). InSAR measurement of the deformation around Siling Co Lake: Inferences on the lower crust viscosity in central Tibet. *Journal of Geophysical Research: Solid Earth*, 120(7), 5290–5310. https://doi.org/10.1002/2014JB011768

Duvall, A. R., & Clark, M. K. (2010). Dissipation of fast strike-slip faulting within and beyond northeastern Tibet. *Geology*, 38(3), 223–226. https://doi.org/10.1130/G30711.1

Dziewonski, A. M., Chou, T.-A., & Woodhouse, J. H. (1981). Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. *Journal of Geophysical Research*, 86(B4), 2825–2852. https://doi.org/10.1029/JB086iB04p02825

Ekström, G., Nettles, M., & Dziewoński, A. (2012). The global CMT project 2004–2010: Centroid-moment tensors for 13,017 earthquakes. *Physics of the Earth and Planetary Interiors*, 200, 1–9. https://doi.org/10.1016/j.pepi.2012.04.002

Elliott, J. R., Walters, R. J., & Wright, T. J. (2016). The role of space-based observation in understanding and responding to active tectonics and earthquakes. *Nature Communications*, 7, 1–16. https://doi.org/10.1038/ncomms13844

England, P., & Houseman, G. (1985). Role of lithospheric strength heterogeneities in the tectonics of Tibet and neighbouring regions. *Nature*, 315(6017), 297–301. https://doi.org/10.1038/315297a0

England, P., & Molnar, P. (1997). Active deformation of Asia: From kinematics to dynamics. *Science*, 278(5338), 647–650. https://doi.org/10.1126/science.278.5338.647
Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., et al. (2007). The Shuttle radar topography mission. *Reviews of Geophysics*, 45(2), 1–33. https://doi.org/10.1029/2005RG000183

Fattahi, H., & Amelung, F. (2014). InSAR uncertainty due to orbital errors. *Geophysical Journal International*, 199(1), 549–560. https://doi.org/10.1093/gji/ggu276

Flesch, L. M., Haines, A. J., & Holt, W. E. (2001). Dynamics of the India-Eurasia collision zone. *Journal of Geophysical Research*, 106(B8), 16435–16460. https://doi.org/10.1029/2001JB000208

Gan, W., Zhang, P., Shen, Z. K., Niu, Z., Wang, M., Wan, Y., et al. (2007). Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements. *Journal of Geophysical Research*, 112(8), 1–14. https://doi.org/10.1029/2005JB004120

Gao, T., Zhang, Y., Kang, S., Abbott, B. W., Wang, X., Zhang, T., et al. (2021). Accelerating permafrost collapse on the eastern Tibetan Plateau. *Environmental Research Letters*, 16(5), 054023. https://doi.org/10.1088/1748-9326/abf7f0

Gaudemer, Y., Tapponnier, P., Meyer, B., Peltzer, G., Shunmin, G., Zhitai, C., et al. (1995). Partitioning of crustal slip between linked, active faults in the eastern Qilian Shan, and evidence for a major seismic gap, the “Tianzhu gap,” on the Western Haiyuan Fault, Gansu (China). *Geophysical Journal International*, 120(3), 599–645. https://doi.org/10.1111/j.1365-246X.1995.tb01842.x

Gomba, G., Rodriguez Gonzalez, F., & De Zan, F. (2017). Ionospheric phase screen compensation for the Sentinel-1 TOPS and ALOS-2 Scan-SAR modes. *IEEE Transactions on Geoscience and Remote Sensing*, 55(1), 223–235. https://doi.org/10.1109/TGRS.2016.2604461

Gruber, S. (2012). Derivation and analysis of a high-resolution estimate of global permafrost zonation. *The Cryosphere*, 6(1), 221–233. https://doi.org/10.5194/tc-6-221-2012

Haines, A. J., Dimitrova, L. L., Wallace, L. M., & Williams, C. A. (2015). Introduction to the vertical derivatives of horizontal stress (VDoHS) rates., 9–18. https://doi.org/10.1007/978-3-319-21578-5_2

Haines, A. J., & Holt, W. E. (1993). A procedure for obtaining the complete horizontal motions within zones of distributed deformation from the inversion of strain rate data. *Journal of Geophysical Research*, 98(B7), 12057–12082. https://doi.org/10.1029/93jb00892

Haines, A. J., & Wallace, L. M. (2020). New Zealand-wide geodetic strain rates using a physics-based approach. *Geophysical Research Letters*, 47(1). https://doi.org/10.1029/2019GL084606

Han, L., Liu-Zeng, J., Yao, W., Shao, Y., Yuan, Z., & Wang, Y. (2021). Coseismic slip gradient at the Western terminus of the 1920 Haiyuan Mw 7.9 earthquake. *Journal of Structural Geology*, 152, 104442. https://doi.org/10.1016/j.jsg.2021.104442

He, J., Ma, J., Zhang, P., Tian, L., Zhu, G., Mike Edmunds, W., & Zhang, Q. (2012). Groundwater recharge environments and hydrogeochemical evolution in the Jiuquan Basin, Northwest China. *Applied Geochemistry*, 27(4), 866–878. https://doi.org/10.1016/j.apgeochem.2012.01.014

Hussain, E., Wright, T. J., Walters, R. J., Bekaert, D. P. S., Lloyd, R., & Hooper, A. (2018). Constant strain accumulation rate between major. *Nature Communications*, 9(1392), 1–9. https://doi.org/10.1038/s41467-018-03739-2

Institute of Geology, China Earthquake Administration & Ningxia Bureau of China Earthquake Administration (IGCEA & NBCEA). (1990). *Active Haiyuan fault zone monograph*. Seismological Press.

International Seismological Centre (ISC). (2019). *ISC-GEM Earthquake Catalogue*. https://doi.org/10.31905/D808B825

Jolivet, R., Candela, T., Lasserre, C., Renard, F., Klinger, Y., & Doin, M. P. (2015). The burst-like behavior of aseismic slip on a rough fault: The creeping section of the Haiyuan fault, China. *Bulletin of the Seismological Society of America*, 105(1), 480–488. https://doi.org/10.1785/0120140237

Jolivet, R., Lasserre, C., Doin, M. P., Peltzer, G., Avouac, J. P., Sun, J., & Dailu, R. (2013). Spatio-temporal evolution of aseismic slip along the Haiyuan fault, China: Implications for fault frictional properties. *Earth and Planetary Science Letters*, 377, 23–33. https://doi.org/10.1016/j.eps1.2013.07.020

Jolivet, R., Lasserre, C., Doin, M.-P., Guilloso, S., Peltzer, G., Dailu, R., et al. (2012). Shallow creep on the Haiyuan fault (Gansu, China) revealed by SAR interferometry. *Journal of Geophysical Research*, 117(B6). https://doi.org/10.1029/2011JB008732

Karra, K., Kontgis, C., Statman-Weil, Z., Mazzariello, J., Mathis, M., & Brumby, S. (2021). Global land use/land cover with Sentinel-2 and deep learning. In *IGARSS 2021-2021 IEEE international geoscience and remote sensing symposium*. IEEE.

King, G., & Nábělek, J. (1985). Role of fault bends in the initiation and termination of earthquake rupture. *Science*, 228(4702), 984–987. https://doi.org/10.1126/science.228.4702.984

Kirby, E., Harkins, N., Wang, E., Shi, X., Fan, C., & Burbank, D. (2007). Slip rate gradients along the eastern Kunlun fault. *Tectonics*, 26(2), 1–16. https://doi.org/10.1029/2006TC002033

Kreemer, C., Blewitt, G., & Klein, E. C. (2014). A geodetic plate motion and global strain rate model. *Geochemistry, Geophysics, Geosystems*, 15(10), 3849–3889. https://doi.org/10.1002/2014GC005407

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

Lehner, B., & Grill, G. (2013). Global river hydrology and network routing: Baseline data and new approaches to study the world's large river systems. *Hydrological Processes*, 27(15), 2171–2186. https://doi.org/10.1002/hyp.9740. (Data is available at www.hydrosheds.org)

Li, B., Yu, Z., Liang, Z., Song, K., Li, H., Wang, Y., et al. (2014). Effects of climate variations and human activities on runoff in the Zoige alpine wetland in the eastern edge of the Tibetan plateau. *Journal of Hydrologic Engineering*, 19(5), 1026–1035. https://doi.org/10.1061/(asce)he.1943-5584.0000868

Li, C., Zhang, P.-z., Yin, J., & Min, W. (2009). Late Quaternary left-lateral slip rate of the Haiyuan fault, northeastern margin of the Tibetan Plateau. *Tectonics*, 28(5). https://doi.org/10.1029/2008TC002302

Li, H., Zhang, Y., Dong, S., Zhang, J., Sun, Y., & Wang, Q. (2020). Neotectonics of the Bailongjiang and Hanan faults: New insights into late Cenozoic deformation along the eastern margin of the Tibetan Plateau. *GSA Bulletin*, 132(9–10), 1–18. https://doi.org/10.1130/b35374.1

Li, J., Zhang, Y., Li, H., Xiong, J., Li, J., & Wu, T. (2016). Revisiting late quaternary slip-rate along the Maqu segment of the eastern Kunlun fault, northeast Tibet. *Acta Seismologica Sinica (English edition)*, 90(2), 486–502. https://doi.org/10.1111/1755-6724.12685

Li, X., Pierce, I. K. D., Bormann, J. M., Hammond, W. C., Zhang, Z., Li, C., et al. (2021). Tectonic deformation of the northeastern Tibetan Plateau and its surroundings revealed with GPS block modeling. *Journal of Geophysical Research: Solid Earth*, 126(5), 1–15. https://doi.org/10.1029/2020JB020733

Li, Y., Jiang, W., Zhang, J., & Luo, Y. (2016). Space geodetic observations and modeling of 2016 Mw 5.9 Menyuan earthquake: Implications on seismogenic tectonic motion. *Remote Sensing*, 8(6), 1–11. https://doi.org/10.3390/rs8060519

Li, Y., Nocquet, J. M., Shan, X., & Song, X. (2021). Geodetic observations of shallow creep on the Laohushan-Haiyuan fault, northeastern Tibet. *Journal of Geophysical Research: Solid Earth*, 126(6), 1–18. https://doi.org/10.1029/2020JB021576

Li, Y., Shan, X., Qu, C., & Wang, Z. (2016c). Fault locking and slip rate deficit of the Haiyuan-Liupanshan fault zone in the northeastern margin of the Tibetan Plateau. *Journal of Geodynamics*, 102, 47–57. https://doi.org/10.1016/j.jog.2016.07.005
Liang, S., Gan, W., Shen, C., Xiao, G., Liu, J., Chen, W., et al. (2013). Three-dimensional velocity field of present-day crustal motion of the Tibetan Plateau derived from GPS measurements. *Journal of Geophysical Research: Solid Earth*, 118(10), 5722–5732. https://doi.org/10.1002/2013JB010503

Liu-zeng, J., Shao, Y., Klinger, Y., Xie, K., Yuan, D., & Lei, Z. (2015). Variability in magnitude of paleoearthquakes revealed by trenching and historical records, along the Haiyuan Fault, China. *Journal of Geophysical Research: Solid Earth*, 120(12), 8304–8333. https://doi.org/10.1002/2015JB012163

López-Quiroz, P., Doin, M. P., Tupin, F., Briole, P., & Nicolas, J. M. (2009). Time series analysis of Mexico City subsidence constrained by radar interferometry. *Journal of Applied Geophysics*, 69(1), 1–15. https://doi.org/10.1016/j.jappgeo.2009.02.006

Loveless, J. P., & Meade, B. J. (2010). Geodetic imaging of plate motions, slip rates, and partitioning of deformation in Japan. *Journal of Geophysical Research*, 115(B2), 1–35. https://doi.org/10.1029/2008jb006248

Matrau, R., Klinger, Y., Woerd, J. V. D., Li, Z., Xu, X., & Zheng, R. (2019). Late Pleistocene—Holocene slip rate along the Hasi Shan restraining bend of the Haiyuan Fault: Implication for faulting dynamics of a complex fault system. *Tectonics*, 38(12), 4127–4154. https://doi.org/10.1029/2019TC005488

Maubant, L., Pathier, E., Daout, S., Radiguet, M., Doin, M.-P., Kazachkina, E., et al. (2020). Independent component analysis and parametric approach for source separation in InSAR time series at regional scale: Application to the 2017–2018 slow slip event in Guerrero (Mexico). *Journal of Geophysical Research: Solid Earth*, 125(3), e2019JB018187. https://doi.org/10.1029/2019JB018187

McCaffrey, R., King, R. W., Payne, S. J., & Lancaster, M. (2013). Active tectonics of northwestern US inferred from GPS-derived surface velocities. *Journal of Geophysical Research: Solid Earth*, 118(2), 709–723. https://doi.org/10.1029/2012JB009473

Meade, B. J., & Hager, B. H. (2005). Block models of crustal motion in southern California constrained by GPS measurements. *Journal of Geophysical Research*, 110(3), 1–19. https://doi.org/10.1029/2004JB003209

Molnar, P., & Tapponnier, P. (1975). Cenozoic tectonics of Asia: Effects of a continental collision. *Science*, 189(4201), 419–426. https://doi.org/10.1126/science.189.4201.419

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), 5–8. https://doi.org/10.3390/rs12030424

Murphy, B., Müller, S., & Yurchak, R. (2021). GeoStat-framework/PyKrige: v1.6.1 (v1.6.1). Zenodo. https://doi.org/10.5281/zenodo.5380342

Ou, Q., Kulikova, G., Yu, J., Elliott, A., Parsons, B., & Walker, R. (2020). Magnitude of the 1920 Haiyuan earthquake reestimated using seismological and geomorphological methods. *Journal of Geophysical Research: Solid Earth*, 125(8). https://doi.org/10.1029/2019jb019244

Pagani, C., Bodin, T., Métois, M., & Lasserre, C. (2021). Bayesian estimation of surface strain rates from global navigation satellite system measurements: Application to the southwestern United States. *Journal of Geophysical Research: Solid Earth*, 126(6), 1–25. https://doi.org/10.1029/2021JB021905

Pichon, X. L., Fournier, M., & Jolivet, L. (1992). Kinematics, topography, shortening, and extrusion in the India-Eurasia collision. *Tectonics*, 11(6), 1085–1098. https://doi.org/10.1029/92tc01566

Qiao, X., Qu, C., Shan, X., Zhao, D., & Liu, L. (2021). Interseismic slip and coupling along the Haiyuan fault zone constrained by InSAR and GPS measurements. *Remote Sensing*, 13(16), 1–20. https://doi.org/10.3390/rs13163333

Quigley, M. C., Jiménez, A., Duffy, B., & King, T. R. (2019). Physical and statistical behavior of multifault earthquakes: Darfield earthquake case study, New Zealand. *Journal of Geophysical Research: Solid Earth*, 124(5), 4788–4810. https://doi.org/10.1029/2019JB017508

Ren, Z., Zhang, Z., Chen, T., Yan, S., Yin, J., Zhang, P., et al. (2016). Clustering of offsets on the Haiyuan fault and their relationship to paleoearthquakes. *Bulletin of the Geological Society of America*, 128(1–2), 3–18. https://doi.org/10.1130/B31155.1

Rousset, B., Jolivet, R., Simons, M., Lasserre, C., Riel, B., Milillo, P., et al. (2016). An aseismic slip transient on the North Anatolian Fault. *Geophysical Research Letters*, 43(7), 3254–3262. https://doi.org/10.1002/2016GL068250

Sandwell, D. T., & Wessel, P. (2016). Interpolation of 2-D vector data using constraints from elasticity. *Geophysical Research Letters*, 43(10), 10703–10709. https://doi.org/10.1002/2016GL070340

Sathiakumar, S., & Barbot, S. (2021). The stop-start control of seismicity by fault bends along the Main Himalayan Thrust. *Communications Earth & Environment*, 2(1), 1–11. https://doi.org/10.1038/s43247-021-00153-3

Savage, J. C., & Burford, R. O. (1973). Geodetic determination of relative plate motion in central California. *Journal of Geophysical Research*, 78(5), 832–845. https://doi.org/10.1029/jb078i005p00832

Savage, J. C., Gan, W., & Svarc, J. L. (2001). Strain accumulation and rotation in the eastern California shear zone. *Journal of Geophysical Research*, 106(B10), 21995–22007. https://doi.org/10.1029/2000jb000127

Schwartz, D. P., & Sibson, R. H. (1989). Fault segmentation and controls of rupture initiation and termination. *US Geological Survey Open File Report No. 89315*, 458

Shao, Y., Liu-zeng, J., Woerd, J. V. D., Klinger, Y., Oskin, M. E., Zhang, J., et al. (2020). Late Pleistocene slip rate of the central Haiyuan fault constrained from optically stimulated luminescence, 14 C, and cosmogenic isotope dating and high-resolution topography. *GSA Bulletin*, 133(7–8), 1347–1369. https://doi.org/10.1130/B35571.1

Shen, L. (2020). Characterising seismic hazard with InSAR measurements: Cases over large length scales (Unpublished doctoral dissertation). PhD thesis. University of Leeds.

Shen, L., Hooper, A., & Elliott, J. (2019). A spatially varying scaling method for InSAR tropospheric corrections using a high-resolution weather model. *Journal of Geophysical Research: Solid Earth*, 124(4), 1–18. https://doi.org/10.1029/2018JB016189

Shen, L., Hooper, A., Elliott, J., & Wright, T. (2021). Large-scale interseismic deformation along the Altyn Tagh fault determined from Sentinel-1 InSAR. In *Egu general assembly*, (pp. EGU21–10278).

Shen, Z. K., Jackson, D. D., & Ge, B. X. (1996). Crustal deformation across and beyond the Los Angeles basin from geodetic measurements. *Journal of Geophysical Research: Solid Earth*, 101(12), 27957–27980. https://doi.org/10.1029/96jb02544

Smith, W. H., & Wessel, P. (1990). Gridding with continuous curvature splines in tension. *Geophysics*, 55(3), 293–305. https://doi.org/10.1190/1.1442837

Song, X., Jiang, Y., Shan, X., Gong, W., & Qu, C. (2019). A fine velocity and strain rate field of present-day crustal motion of the northeastern Tibetan Plateau inverted jointly by InSAR and GPS. *Remote Sensing*, 11(4), 435. https://doi.org/10.3390/rs11040435

Sudhaus, H., & Jónsson, S. (2009). Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors: Application to the June 2000 Kleifarvatn earthquake, Iceland. *Geophysical Journal International*, 176(2), 389–404. https://doi.org/10.1111/j.1365-246X.2008.03989.x

Tapponnier, P., Zhiqin, X., Roger, F., Meyer, B., Arnaud, N., Wittlinger, G., & Jingsui, Y. (2001). Oblique stepwise rise and growth of the Tibet Plateau. *Science*, 294(5547), 1671–1677. https://doi.org/10.1126/science.105978

Tarantola, A. (2005). Inverse problem theory and methods for model parameter estimation. *Society for Industrial and Applied Mathematics*. https://doi.org/10.1137/1.9780898717921
Thollard, F., Clesse, D., Doin, M.-p., Donadieu, J., Durand, P., Grandin, R., et al. (2021). Flatsim: The ForM@Ter LArge-scale multi-temporal Sentinel-1 InterferoMetric Service. *Remote Sensing*, 13(18), 1–29. https://doi.org/10.3390/rs13183734

Turcotte, D., & Schubert, G. (2014). *Geodynamics* (3rd edn). Cambridge University Press. https://doi.org/10.1017/CBO9780511843877

Tymofeyeva, E., & Fialko, Y. (2018). Geodetic evidence for a blind fault segment at the southern end of the San Jacinto Fault Zone. *Journal of Geophysical Research: Solid Earth*, 123(1), 878–891. https://doi.org/10.1002/2017JB014477

Walters, R. J., Gregory, L. C., Wedmore, L. N., Craig, T. J., McCaffrey, K., Wilkinson, M., et al. (2018). Dual control of fault intersections on stop-start rupture in the 2016 Central Italy seismic sequence. *Earth and Planetary Science Letters*, 500, 1–14. https://doi.org/10.1016/j.eps1.2018.07.043

Wang, H., Liu-Zeng, J., Ng, A.-M., Ge, L., Javed, F., Long, F., et al. (2017). Sentinel-1 observations of the 2016 Menyuan earthquake: A buried reverse event linked to the left-lateral Haiyuan fault. *International Journal of Applied Earth Observation and Geoinformation*, 61, 14–21. https://doi.org/10.1016/j.jag.2017.04.011

Wang, H., & Wright, T. J. (2012). Satellite geodetic imaging reveals internal deformation of western Tibet. *Geophysical Research Letters*, 39(7), 1–5. https://doi.org/10.1029/2012GL051222

Wang, H., Wright, T. J., Liu-Zeng, J., & Peng, L. (2019). Strain rate distribution in south-central Tibet from two decades of InSAR and GPS. *Geophysical Research Letters*, 46(10), 5170–5179. https://doi.org/10.1029/2019GL081916

Wang, J., Wang, C., Zhang, H., Tang, Y., Duan, W., & Dong, L. (2021). Freeze-thaw deformation cycles and temporal-spatial distribution of permafrost along the Qinghai-Tibet railway using multitrack InSAR processing. *Remote Sensing*, 13(23), 4744. https://doi.org/10.3390/rs13234744

Wang, M., & Shen, Z.-K. (2020). Present-day crustal deformation of continental China derived from GPS and its tectonic implications. *Journal of Geophysical Research: Solid Earth*, 125(2), 1088–1102. https://doi.org/10.1093/gji/ggw445

Wang, Y., Chang, L., Feng, W., Samsonov, S., & Zheng, W. (2021). Topography-correlated atmospheric signal mitigation for InSAR applications in the Tibetan plateau based on global atmospheric models. *International Journal of Remote Sensing*, 42(11), 4361–4379. https://doi.org/10.1080/01431161.2021.1892856

Wang, Y., Wang, M., & Shen, Z. K. (2017). Block-like versus distributed crustal deformation around the northeastern Tibetan plateau. *Journal of Asian Earth Sciences*, 140, 31–47. https://doi.org/10.1016/j.jseas.2017.02.040

Weiss, J. R., Walters, R. J., Morishita, Y., Wright, T. J., Lazecky, M., Wang, H., et al. (2020). High-resolution surface velocities and strain for Anatolia from Sentinel-1 InSAR and GNSS data. *Geophysical Research Letters*, 47(17). https://doi.org/10.1029/2020GL087376

Wesnousky, S. G. (2008). Displacement and geometrical characteristics of earthquake surface ruptures: Issues and implications for seismic-hazard analysis and the process of earthquake rupture. *Bulletin of the Seismological Society of America*, 98(4), 1609–1632. https://doi.org/10.1785/0120070111

Wessel, P., Luis, J. F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W. H., & Tian, D. (2019). The generic mapping tools version 6. *Geochemistry, Geophysics, Geosystems*, 20(11), 5556–5564. https://doi.org/10.1029/2019GC008515

Wessel, P., Smith, W. H. F., Scharroo, R., Luis, J., & Wobbe, F. (2013). Generic mapping tools: Improved version released. *Eos, Transactions American Geophysical Union*, 94(45), 409–410. https://doi.org/10.1002/2013EO450001

Wicks, C., Weaver, C., Bodin, P., & Sherrod, B. (2013). InSAR Evidence for an active shallow thrust fault beneath the city of Spokane Washington, USA. *Journal of Geophysical Research: Solid Earth*, 118(3), 1268–1276. https://doi.org/10.1002/jgrb.50118

Wright, T. J., Elliott, J. R., Wang, H., & Ryder, I. (2013). Earthquake cycle deformation and the Moho: Implications for the rheology of continental lithosphere. *Tectonophysics*, 609, 504–523. https://doi.org/10.1016/j.tecto.2013.07.029

Wu, H., Wang, L., Liu, L., Hao, L., Ma, Q., Li, F., et al. (2019). Groundwater recharge and evolution in the Wuwei Basin, northwestern China. *Environmental Earth Sciences*, 78(12), 1–10. https://doi.org/10.1007/s12665-019-8362-5

Xiong, X., Shan, B., Zheng, Y., & Wang, R. (2010). Tectonophysics stress transfer and its implication for earthquake hazard on the Kunlun fault, Tibet. *Tectonophysics*, 482(1–4), 216–225. https://doi.org/10.1016/j.tecto.2009.07.020

Xu, X., Han, Z., Yang, X., Zhang, S., Yu, G., Zhou, B., et al. (2016). *Zhong Guo Ji Lin Jin Di Qu Di Zhen Gou Zao Tu*. https://doi.org/10.12031/activefault.china.250.2016

Xu, X., & Sandwell, D. T. (2020). Toward absolute phase change recovery with InSAR: Correcting for Earth tides and phase unwrapping ambiguities. *IEEE Transactions on Geoscience and Remote Sensing*, 58(1), 726–733. https://doi.org/10.1109/tgrs.2019.2940207

Xu, X., Sandwell, D. T., Klein, E., & Bock, Y. (2021). Integrated sentinel-1 InSAR and GNSS time-series along the San Andreas fault system. *Journal of Geophysical Research: Solid Earth*, 126(11), 1–14. https://doi.org/10.1029/2021JB022579

Xu, X., Zhang, Z., Hu, F., & Chen, X. (2019). Dynamic rupture simulations of the 1920 Ms 8.5 Haiyuan earthquake in China. *Bulletin of the Seismological Society of America*, 109(5), 2009–2020. https://doi.org/10.1785/0120190061

Yao, W., Liu-Zeng, J., Oskin, M. E., Wang, W., Li, Z., Prush, V., et al. (2019). Reevaluation of the late Pleistocene slip rate of the Haiyuan fault near Songshan, Gansu province, China. *Journal of Geophysical Research: Solid Earth*, 124(5), 5217–5240. https://doi.org/10.1029/2018JB016907

Yu, C., Li, Z., & Blewitt, G. (2021). Global comparisons of ERA5 and the operational HRES tropospheric delay and water vapor products with GPS and MODIS. *Earth and Space Science*, 8(5), 1–18. https://doi.org/10.1029/2020EA001417

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

Yuan, D. Y., Ge, W. P., Chen, Z. W., Li, C. Y., Wang, Z. C., Zhang, H. P., et al. (2013). The growth of northeastern Tibet and its relevance to large-scale continental geodynamics: A review of recent studies. *Tectonics*, 32(5), 1358–1370. https://doi.org/10.1002/tect.20081

Zhang, J., Bock, Y., Johnson, H., Fang, P., Williams, S., Genrich, J., et al. (1997). Southern California Permanent GPS Geodetic Array: Error analysis of daily position estimates and site velocities. *Journal of Geophysical Research*, 102(B8), 18035–18005. https://doi.org/10.1029/97JB01380

Zhang, J., Ma, Z., Li, Z., & Li, W. (2009). Study on deformational characteristics of northeastern Qinghai-Xizang (Tibetan) Plateau from late cenozoic deformation in the Xining Basin. *Geological Review*, 55(4), 731–747.

Zhang, P., Molnar, P., Burchfiel, B. C. C., Royden, L., Wang, Y., Deng, Q., et al. (1988). Bounds on the Holocene slip rate of the Haiyuan fault, north-central China. *Quaternary Research*, 30(2), 151–164. https://doi.org/10.1016/0033-5894(88)90020-8

Zhang, X., Wang, N., Xie, Z., Ma, X., & Huete, A. (2018). Water loss due to increasing planted vegetation over the Badain Jaran Desert, China. *Remote Sensing*, 10(1), 1–21. https://doi.org/10.3390/rs10010134

Zhang, Y., Shan, X., Zhang, G., Zhong, M., Zhao, Y., Wen, S., et al. (2020). The 2016 Mw 5.9 Menyuan earthquake in the Qilian Orogen, China: A potentially delayed depth-segmented rupture following from the 1986 Mw 6.0 Menyuan earthquake. *Seismological Research Letters*, 91(2), 758–769. https://doi.org/10.1785/0220190168
Zhao, S., Zeng, R., Zhang, H., Meng, X., Zhang, Z., Meng, X., et al. (2022). Impact of water level fluctuations on landslide deformation at Long-
yangxia reservoir, Qinghai province, China. *Remote Sensing, 14*(212), 1–33. https://doi.org/10.3390/trs14010212

Zheng, G., Wang, H., Wright, T. J., Lou, Y., Zhang, R., Zhang, W., et al. (2017). Crustal deformation in the India-Eurasia collision zone from 25 Years of GPS measurements. *Journal of Geophysical Research: Solid Earth, 122*(11), 9290–9312. https://doi.org/10.1002/2017JB014465

Zheng, W., Yuan, D., Zhang, P., Yu, J., Lei, Q., Wang, W., et al. (2016). Tectonic geometry and kinematic dissipation of the active faults in the northeastern Tibetan Plateau and their implications for understanding northeastern growth of the plateau. *Quaternary Sciences, 36*(4), 775–788. https://doi.org/10.11928/j.issn.1001-7410.2016.04.01

Zheng, W., Zhang, P., He, W., Yuan, D., Shao, Y., Zheng, D., et al. (2013). Transformation of displacement between strike-slip and crustal shortening in the northern margin of the Tibetan Plateau: Evidence from decadal GPS measurements and late Quaternary slip rates on faults. *Tectonophysics, 584*, 267–280. https://doi.org/10.1016/j.tecto.2012.01.006

Zheng, W. J., Zhang, P. Z., Ge, W. P., Molnar, P., Zhang, H. P., Yuan, D. Y., & Liu, J. H. (2013). Late quaternary slip rate of the South Heli Shan fault (northern Hexi corridor, NW China) and its implications for northeastward growth of the Tibetan plateau. *Tectonics, 32*(2), 271–293.
https://doi.org/10.1002/tect.20022

Zhou, Y., Thomas, M. Y., Parsons, B., & Walker, R. T. (2018). Time-dependent postseismic slip following the 1978 Mw 7.3 Tabas-e-Golshan, Iran earthquake revealed by over 20 years of ESA InSAR observations. *Earth and Planetary Science Letters, 483*, 64–75. https://doi.org/10.1016/j.epsl.2017.12.005

Zhu, Y., Wang, K., & He, J. (2020). Effects of earthquake recurrence on localization of interseismic deformation around locked strike-slip faults. *Journal of Geophysical Research: Solid Earth, 125*(8). https://doi.org/10.1029/2020JB019817

Zou, D., Zhao, L., Sheng, Y., Chen, J., Hu, G., Wu, T., et al. (2017). A new map of permafrost distribution on the Tibetan Plateau. *The Cryosphere*, *11*(6), 2527–2542. https://doi.org/10.5194/tc-11-2527-2017