Asymmetric Interseismic Strain across the Western Altyn Tagh Fault from InSAR

Abstract

As the northern boundary of the Tibetan Plateau, the long Altyn Tagh fault (ATF) controls the regional tectonic environment, and the study of its long-term fault slip rate is key to understanding the tectonic evolution and deformation of the northern Tibetan Plateau. In this paper, we measure the fault slip rate of the western segment of the ATF using InSAR observations between 2015 to 2020. The Multi-Temporal Interferometric InSAR analysis is applied to obtain the two-dimensional fault-parallel and vertical displacement fields. The spatially dense InSAR observations clearly illustrate the asymmetrical pattern of displacement fields across the fault. Constrained by our InSAR observations, the fault slip rate and locking depth of the western segment of the ATF are inverted using four different models in a Bayesian framework. The two-layer viscoelastic model incorporating lateral heterogeneity of rheology in the lower crust indicates that the fault slip rate of the western ATF is estimated to be 9.8 ± 1.1 mm/yr (at 83.8°E across the ATF) and 8.6 ± 1.1 mm/yr (at 85.1°E), respectively, and the locking depth is 15.8 ± 4.3 km and 14.8 ± 4.9 km. Our new estimates generally agree with the previous estimates of fault slip rate constrained by GPS observations. We conclude that the contrast between the thickness of the elastic layer and the shear modulus of the Tibetan plateau and the Tarim basin jointly contribute to the asymmetric interseismic strain accumulation on the ATF.

1. Introduction

The 1600-km-long Altyn Tagh fault (ATF) is the northern boundary of the Tibetan Plateau (TP). It separates the rigid Tarim basin from the relatively weak TP [1,2] and connects almost all orogenic belts and the thrust systems on the northern margin of the plateau. Investigating the earthquake cycle behavior of this large left-lateral strike-slip fault system informs us of the uplift process, stress transfer, and crustal deformation mechanism of the TP [3,4,5]. For instance, the slip rate of the ATF has been a matter of debate to discriminate among continental deformation models [6,7,8]. The large fault slip rate (20–30 mm/yr) of the ATF supports the ‘lithosphere extrusion model,’ which proposes that the India–Eurasia collision be accommodated laterally along large-scale narrow bounding fault zones by lithosphere extrusion [2,9,10]. Whereas slow slip rate along the ATF supports the ‘continuum model,’ which stresses that the India–Eurasia collision is partitioned by numerous, isolated fault structures [9,10,11]. Determining the fault slip rate of ATF helps to understand the interrelationship between the large-scale plateau boundary faults and continental dynamics [12,13,14].

Due to the harsh environment and inconvenient transportation in most areas along with the ATF, conventional geodetic measurements such as Global Positioning System (GPS) and leveling are limited. Interferometric Synthetic Aperture Radar (InSAR) substantially improves the temporal and spatial resolution of surface deformation via stacking individual interferograms, thus providing crustal deformation maps at a decadal scale [15,16,17,18]. Wright et al. [19] processed the European Remote Satellite (ERS) images (1992–1999) across the western ATF (79°E) and obtained the interseismic fault slip rates of the Karakoram fault and the western ATF. Their work pioneered the use of InSAR to obtain interseismic crustal deformation. Elliott et al. [20] found that the ATF (85°E) slips at 11 ± 5 (1σ) mm/yr, and their derived line-of-sight (LOS) displacement field features an asymmetric pattern across the fault. However, due to the phase decorrelation in interferograms around the Altyn mountains, the LOS velocities show a sense of asymmetry to the north of the ATF, opposite to that expected for the rigid Tarim basin. Zhu et al. [21] processed the Environmental Satellite (ENVISAT) data spanning 2003–2010 at almost the same location as Elliott et al. [20] but identified an asymmetric velocity pattern to the south of the ATF. Jolivet et al. [22] processed InSAR images across the eastern ATF (94°E) and observed apparent velocity asymmetry on both sides of the fault. They linked the asymmetric strain pattern to rigidity differences across the ATF or the geometry of a south-dipping fault. More recently, using a combination of GPS and InSAR measurements at 86°E, Li et al. [23] favored no asymmetric strain pattern across the ATF. Several factors might contribute to the discrepancies among previous studies. For example, vertical crustal deformation is neglected in the InSAR observations, which are more sensitive to vertical deformation than horizontal deformation. In addition, simplified two-dimensional elastic models might not be reasonable for modeling the interseismic strain.

Compared with previous ERS and ENVISAT satellites, the Sentinel-1 satellite has the advantages such as wide coverage, a short revisit period, and free public access to the data policy. Besides, Terrain Observation with Progressive Scans (TOPs) imaging technology is adopted, ensuring a large area of ground monitoring and reducing the revisit period [24,25,26]. The Sentinel SAR data have been widely used in geological disaster monitoring [27,28,29].

In this study, we processed three ascending and three descending Interferometric Wide swath (IW) mode images across the western ATF. We decomposed the ascending and descending LOS maps to obtain the horizontal and vertical deformation fields. We modeled the fault-parallel velocity profiles using different models to investigate possible asymmetric interseismic strain across the ATF. Our results provide new insights into the kinematic characters along with the ATF.

2. Data and Methods

2.1. GPS Data

He et al. [30] established a 400-km-long transect consisting of 17 GPS stations across the ATF (86°E) in 2009 (indicated by the red triangle symbol in Figure 1). Those campaign GPS stations was routinely observed for 2–3 periods in 2009–2011. In December 2017, Li et al. [31] re-measured the same GPS stations and processed those campaign-mode GPS data together with data at seven regional Crustal Movement Observation Network of China (CMONOC, https://iag.dgfi.tum.de/media/archives/prchina03/CRUSTAL%20MOVEMENT%20OBSERVATION.htm (accessed on 8 February 2022)) continuous sites near this region. We combined the GPS velocity field of Li et al. [31] with that of Wang and Shen [32] by minimizing the velocity residuals within common sites. The combined GPS velocity is referenced to a Eurasian frame. To highlight the relative motion of the crust on both sides of the Altyn Tagh fault, we transformed the GPS site velocities into a Tarim-fixed reference frame. The Euler pole in Wang et al. [33] is selected, then the GPS velocity field relative to Tarim Basin is obtained by deducting the rotational motion and global motion of the Euler pole from the whole velocity field data. A more detailed description of the GPS data process can be found in Li et al. [31].

2.2. InSAR Data and Processing

Three ascending and three descending Interferometric Wide (IW) swath mode imagery are selected, and three frames on each orbit are concentrated into a long strip. The coverage of each strip is about 250 km × 700 km, and the observation period is from 2014 to 2020. See Figure 1 for InSAR data coverage.

GAMMA software performs the InSAR processing [35], starting with Single Look Complex (SLC) level data. For each track, we take the middle acquisition as the single prime, then use the imaging geometry as the reference for the images of other dates, and perform registration resampling after stitching until the registration accuracy of all images relative to the reference is higher than 1/1000 pixels in the azimuth direction [36,37]. In this way, interferograms can be generated with all combinations among those scenes without re-registration. The selection criterion of the interferogram pairs according to the spatial baseline is less than 30 m, and the time baseline is more than 1 year and less than 3 years, considering that the study area in the summer by the influence of the thickness of active layers change with the melting permafrost [38,39,40,41,42], to eliminate the interferograms between June to September [43,44], the final baseline network is shown in Figure A1. In the interferometry process, a multilooking operation with 20 looks in the range direction and 4 looks in the azimuth direction is applied to suppress noises. The simulation and removal of the topography phase are performed using the 30 m ALOS Global Digital Surface Model ALOS World 3d-30 m (AW3D30) [45,46], and the two-pass InSAR approach is used to produce differential interferograms, then the interferograms are filtered using a power spectrum filter [47]. The rock area with a stable phase is selected as a reference point and the region with coherence smaller than 0.2 is masked. The unwrapped phase is derived from the filtered phase using a minimum cost flow (MCF) algorithm [48]. Finally, the unwrapped phase is geocoded in the geographic coordinate system.

For ATF, an area with significant relief, spatial variations in hydrostatic delays are strong, it is necessary to reduce the tropospheric effects carefully [49,50,51]. The traditional atmospheric correction method based on the empirical linear correlation formula of phase elevation may remove the concerned deformation signal [52,53,54,55], especially in this region where long-wavelength deformation signals correlate with topography. For this purpose, the stratified tropospheric delay is corrected by Generic Atmospheric Correction Online Service (GACOS) [56], which is based on reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF) with a spatial resolution of 0.125° and a temporal resolution of 6 h. The Iterative Tropospheric Decomposition (ITD) model is used to separate the elevation-related signals and turbulence signals from the Zenith Tropospheric Delay (ZTD) and interpolate this to generate high-resolution tropospheric delay maps, and the high-resolution tropospheric delay is obtained by interpolation [56,57,58]. Compared with other atmospheric delay correction methods, its effectiveness has been verified [59,60,61], and it has been applied in many studies recently and achieved good results [62,63,64].

Detailed atmospheric phase correction steps of GACOS can be found in reference [65]. The interferograms corrected by the GACOS atmospheric correction are obtained by further differencing between the unwrapped interferograms and the atmospheric correction maps. A quadratic surface model and iterative least-squares method are used to remove residual orbit errors. Here, the correction process is briefly described. Figure 2 shows the 12-day ascending track interferogram from 15 to 27 December 2019 (20191215–20191227). Since the interval is only 12 days, it can be considered that the interferogram does not contain a crustal deformation signal, while atmospheric behavior and orbital error play a dominant role. Figure 2b shows the phase difference between two GACOS ZTD maps obtained at the date of Sentinel-1 acquisitions and project the result into line-of-sight (LOS) using the incidence angle of the SAR image. By only applying the GACOS corrections, the corrected phase still shows an insignificant correlation with the topography, we observed long-wavelength residual signals, and we also used the spatial information contained in the cumulative displacement to determine the empirical correction on each interferogram [66]. Figure 2 shows an example where we have removed a quadratic ramp in latitude and a linear ramp in longitude, a latitudinal cross dependency of the phase-elevation relationship, and then apply the procedure to all the interferograms. Figure 2d shows the phase value after correction. From the histogram, the original interference phase diagram and GACOS simulation have two peaks (Figure 2g,h). After correction, the phase histogram has only one peak value and is near 0 value, showing obvious atmospheric delays have been most effectively corrected. However, we admit that there is still a left side heavy tail distribution in the GACOS corrected phases, as shown in Figure 2i. As can be seen from the original interferogram, this signal has always existed in the upper left corner of the graph, and we interpret it as atmospheric turbulence signal, so neither GACOS nor the phase-elevation empirical formula can remove it, but it can be removed by filtering in the subsequent time-series processing.

The corrected interferograms are then imported into the poly-interferogram rate and time-series estimator algorithm (Π-RATE) software package for time-series analysis [67,68,69]. Unlike other Small Baseline Subset (SBAS) methods, Π-RATE uses a minimum distance spanning tree algorithm to pick out independent observations for time series analysis. In this method, pixels that are not always coherent in different interferograms are employed to calculate deformation rates to provide as much constraint on the rate as possible [70,71]. The observed values also include a random Atmospheric Phase Screen (APS), a temporal high-pass filter followed by a spatial low-pass filter applied to mitigate atmospheric disturbances. LOS velocity and its uncertainty are obtained using an iterative weighted least squares approach on a pixel-wise basis.

3. Results

3.1. InSAR LOS Velocity Field

To obtain a large-scale velocity field including several tracks, it is a common procedure that ties the InSAR to the GNSS velocities like Weiss et al. [72]. However, due to the scarity of GPS stations (most of the research area is uninhabited), the final result of tie to GPS may be biased by removing the best-fit plane. Here we adopted a simple method; first, we calculate the datum offset using the overlapped zone of the adjacent track, then the offset is subtracted from the adjacent track to unify the datum. After removing the average offset in each overlapping area between every two adjacent tracks in ascending and descending, respectively, we obtain the merged LOS velocity fields. We project the GPS velocities into the local satellite line-of-sight and calculate the difference from the InSAR velocities. Three GPS stations (AT13, AT14, and AT15) are chosen as reference sites to align those datasets, as shown in Figure 3a,c. all the velocity maps were computed with respect to an average value at the location of the three GPS stations. There is an obvious cross-fault deformation gradient on both sides of the fault, and the color is relatively uniform on the south or north side of the fault, but there is an obvious difference near the fault zone. The blue movement is toward the satellite, and the red movement is away from the satellite, which is consistent with a left-lateral motion across the fault.

Figure 3b,d shows the cross-fault rate profile (PQ) in the ascending and descending InSAR rate maps. From the profile of PQ, there are obvious velocity gradients across the fault, and the LOS velocity difference is about 4–5 mm/yr. In the study area, the ATF is the major fault, whereas the velocity profiles also reflect the deformation of the Altyn North Margin fault (NATF) and Cherchen fault (CHEF). Such an effect may lead to an asymmetrical pattern in the velocity profiles. The ascending and descending velocity fields exhibit a clear transition, reflecting the left-lateral strike-slip fault motion.

We compare the LOS-converted GPS (no vertical velocities) profile (Figure 1) in He et al. [30] at longitude 86°E (Figure 3c,d) with our InSAR observations. The result shows that the difference between them is distributed within ±1 mm/yr, comparable to InSAR measurement accuracy (about 1 mm/yr). The good agreement indicates the reliability of our InSAR rate map.

3.2. ATF-Parallel and Vertical Surface Velocities

According to the previous tectonic results, considering the geometry, structural continuity of ATF, and the earthquakes, the whole Altyn Tagh fault system can be divided into eight segments [73]. We are concerned about the Nierkule and Kurukuole segments in this paper. The continuous Nierkule section from 82.75°E to 84.5°E has a length of 190 km, characterized by the development of a fault depression valley. The Kurukuole segment from 84.5°E to 87.2°E has a length of about 230 km and an overall trend of NE70°, characterized by obvious linear characteristics and a large fault slip rate.

When InSAR observations are made in the northeastern margin of Qinghai-Tibet Plateau, LOS deformation inherently contains a contribution from vertical deformation, which may lead to deviation of crustal deformation interpretation based on InSAR data. Therefore, after obtaining the velocity fields from different perspectives of the ascending and descending orbits, we derive maps of horizontal and vertical average velocities by combining different orbits. By combining two independent observation directions and assuming that the third component is zero, we can project the information in the direction of any two orthogonal basis vectors [74]. Here, we use the method in Lindsey et al. [75]:

$$\left[\begin{array}{cc}{E}_{a}sin\phi +N_{a}cos\phi & U_a\\ E_{d}sin\phi +N_{d}cos\phi & U_d\end{array}\right]\left[\begin{array}{c}{D}_H\\ D_V\end{array}\right]=\left[\begin{array}{c}LO{S}_a\\ LO{S}_d\end{array}\right]$$

(1)

where E, N, and U are east, north, and up components (subscript character a for ascending and d for descending) of the unit look vectors for the satellite tracks considering the radar incidence angle; D_H and D_V are horizontal and vertical surface displacement components; $\phi$ is the local fault strike, and $LO{S}_a$ and $LO{S}_d$ are InSAR LOS displacements from the ascending and descending tracks, respectively. Assuming that the horizontal movement is N70°E along the average strike of ATF, Figure 4 shows the fault-parallel and vertical displacement results.

To compare the decomposed components with GPS results, GPS measurements are projected in parallel and vertical fault directions and compared with InSAR velocity profiles. Therefore, in addition to the PQ profile, an MN profile is added in the west, as shown in Figure 4. Within 100 km on both sides of the MN profile, there are only a few GPS stations on one side of the Tarim Basin.

The most prominent feature in the deformation field is the large-scale tectonic motion across the fault boundary with sharp transitions. Superimposed on this tectonic pattern are numerous nontectonic deformations mostly related to permafrost. Firstly, in addition to the strike-slip rate of the ATF, the Cherchen fault zone is also reflected in the rate map. Cherchen fault zone is a boundary fault in the northwestern part of the southeastern Tarim Basin, which is approximately parallel to the ATF. The fault has been active for a long time since the end of the early Paleozoic and has controlled the sedimentation and affected the tectonic deformation in the southeastern Tarim Basin. The SAR data mainly cover the western segment, with strong Cenozoic activity intensity and developed branch faults, which can be divided into Cherchen deep fault and Cherchen main fault. The western segment of the Cherchen fault bears not only the compressive stress caused by the west Kunlun uplift but also the compressive stress caused by the Altyn uplift, so the fault activity intensity is high, and branch faults are developed. At the same time, because of the multiple compressive stress directions in the western section, the Cherchen fault plane is complex and has the characteristics of strike-slip thrust activity. This characteristic is shown in both vertical and horizontal decomposition rate fields, and similar characteristics are also shown in Daout et al. [40]. In Figure 4b, there is a vertical component, which shows the thrust characteristic.

Secondly, obvious and steep displacement gradient changes are observed on the Manyi fault trace, and the post-earthquake effect can still be observed. The post-earthquake deformation seems to have a local pattern, with the deformation mainly concentrated at about 50 km from the fault. The seismogenic fault of the Manyi earthquake in 1997 is a strike-slip fault with a relatively simple structure. Using InSAR, the surface coseismic rupture length is about 170 km, the maximum displacement is 7 m, and the overall strike is N76°E [76]. Some scholars have studied the post-earthquake mechanism of the Manyi earthquake. Ryder et al. [77] found that both the standard linear body model and the afterslip model can be used to describe the post-earthquake deformation. It has been 24 years since the earthquakes occurred, and further research is needed to determine whether it has gradually entered the interseismic stage from post-earthquake [78,79].

Finally, there are large areas with sporadic permafrost in the vertical rate map. The annual average temperature of the Qinghai–Tibet Plateau is very low, forming a large area of permafrost, which is the highest elevation and largest permafrost area in the world [80,81,82,83]. Permafrost in the Tibetan Plateau is very sensitive to climate change and unstable in nature [84]. In the cold environment, the soil is frozen, the water in the soil is redistributed during the freezing process, and the structure and density of soil particles change with the generation of pressure, resulting in frost heave. In a warm environment, the frozen soil melts, the soil drains and consolidates under the action of gravity and external load, and the soil layer compresses and deforms, resulting in settlement. These are the characteristics of frost heaving and thawing settlement, and such changes will be repeated and reversible with the alternation of seasons [85,86,87]. Thawing subsidence and freezing uplift are reflected in the InSAR rate map [40]. Therefore, when selecting data for time series processing, we try to avoid the summer with serious frozen soil melting to minimize the impact [44]. Nevertheless, there are still large vertical deformation variables near major glaciers and the Yanghu basin.

3.3. Locking Depth and Slip Rate

The screw dislocation model of Savage & Burford [88] in homogeneous elastic half-space is used to retrieve the far-field velocity and the locking depth of ATF firstly. The model assumes that the slip on the fault is zero (locked) to depth D but is an amount $v_0$ everywhere below that depth. The part above locking depth cannot slide freely, the shear strain will accumulate near the fault, and the strain value will increase with the distance closer to the fault. For an infinitely long vertical strike-slip fault, the relationship between slip rate and distance can be expressed by an arctangent function as follows [70]:

$$V_x=\frac{v_0}{\pi }ta{n}^{-1}\left(\frac{x}{D}\right)+V_{offset}$$

(2)

where $V_x$ is observed fault-parallel velocity, $v_0$ is the slip rate, x is the perpendicular distance to the fault trace, D is the locking depth, and $V_{offset}$ is the rate offset constant related to the reference frame [89].

To solve the locking depth and slip rate along with other parameters given in Equation (2) with the uncertainty in consideration, we use the Bayesian probability density method to find the best-fitting fault parameters as well as their uncertainties [89,90,91,92]. Referring to previous research results, a prior constraint is carried out on the range of parameter values. The slip rate is constrained from 1 mm/yr to 30 mm/yr, the locking depth is constrained from 1 km to 65 km, and all parameters are a uniform prior probability distribution over each range. Based on the Markov chain Monte Carlo (MCMC) sampler [93,94], 600 initial walkers are selected to explore the parameter space. For each InSAR slip rate profile across the fault, the model ran over 1 million iterations and produced over 20,000 independent random samples from which we estimate both the maximum a posteriori probability (MAP) solution and corresponding parameter uncertainties. The MAP solution is taken as the optimal solution for the model parameters.

The results of our analysis are shown in Figure 5. The best-fitting results (blue curve in the figure) show that the slip rate of the Nierkule segment is 7 ± 0.1 mm/yr and the locking depth is 18.1 ± 0.6 km. At the Kurukuole segment, the fault slip rate decreases to 6.1 ± 0.3 mm/yr, and the locking depth becomes 11.7 ± 0.4 km. The results show that the slip rate and locking depth change obviously from west to east and decrease gradually from west to east.

4. Discussion

4.1. Strain Asymmetry

Over the last two decades, the proliferation of geodetic measurements, including GPS and InSAR, has helped to document asymmetric patterns of cross-fault velocity on major continental strike-slip faults [95,96,97], such as the San Andreas fault [98,99,100] and the Anatolian fault [101,102]. Elliott et al. [20] and He et al. [30] observed asymmetric patterns in the Altyn fault zone and excluded the possibility of offsetting the location of the fault trace, meaning that the deep fault plane shifted southward for more than 10 km. Therefore, the southward dip of the fault is nearly 45°, which is obviously unreasonable for the strike-slip ATF with a vertical fault plane. Jolivet et al. [22] found that the rigidity contrast coefficient across the fault is 0.85 and concluded that such asymmetry is the joint effect of a rigidity decrease from Tarim to Qaidam and a southward fault trace deviation.

The possible factors leading to the asymmetric velocities across the strike-slip faults are summarized as follows: (1) the difference in medium properties on both sides of the fault [103]; (2) the varying thickness of the elastic layer (hereafter Te) on both sides of the fault [104]; (3) the rheological contrast in the lower crust and upper mantle across faults [105]; (4) the deviation of the fault positions relative to fault trace below the brittle layer [106].

The reason in part accounting for the above disputes is the scarcity of geodetic measurements in some regions and the limited resolution of the diverse models. Constrained by the high-resolution rate map, we establish different models to fit and explain the asymmetry of the velocity profile across the ATF. In Model A, we used the Savage–Burford model in an elastic half-space (Figure 6a); Models B and C have two layers, and they are both the Savage–Prescott coupling model but different in thickness of the elastic layer (Figure 6b,c). Model D is the three-layer channel viscous model (Figure 6d), which is discussed in Section 4.1.3.

4.1.1. Modified Half-Space Model

For the Model A, we calculate an asymmetric surface velocity field on each side of the fault accounting for the medium difference on both sides of the fault [22]:

$$V_x\{\begin{array}{cc}\frac{2K{v}_0}{\pi }{tan}^{-1}\left(\frac{x}{D}\right)+V_{offset}\hfill & \mathrm{if} x>0\hfill \\ \frac{2\left(1-K\right)v_0}{\pi }{tan}^{-1}\left(\frac{x}{D}\right)+V_{offset}\hfill & \mathrm{if} x>0\hfill \end{array}$$

(3)

where K is the rigidity contrast, representing the rigidity ratio of the upper crust on both sides of the fault, ranging from 0 to 1. If K = 0.5, it means that the rigidity of the upper crust on both sides is the same. Inversion results are shown in Figure 7a and Figure 8a. The ratio for rigidity of the upper crust in the Nierkule segment is 0.6/0.4, and in the Kurukuole segment is 0.7/0.3, which means that the crustal rigidity on the north side of the fault is about 0.5–1.3 times larger than that on the south side. At the same time, the resolved slip rate of 6.8 ± 0.1 mm/yr is similar to the inversion result of the model without considering the rigidity parameters. The locking depth of 17.2 ± 0.4 km is about 1 km shallower in the Nierkule segment and 1 km deeper in the Kurukuole segment. Generally, the inversion result with the addition of the asymmetry coefficient does not change much for the elastic half-space model.

4.1.2. Viscoelastic-Coupling Models

Both Model B and C are viscoelastic-coupling models [108]. It is assumed that an elastic upper crust is overlying a viscoelastic (the Maxwell body) lower crust or upper mantle, and there is a linear coupling relationship between the two layers.

$$v\left(x,z\right)=\frac{\Delta u}{\pi t_R}{e}^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{$t_R$}\right.}{\displaystyle \sum }_{n=1}^{\infty }\frac{{\left(\raisebox{1ex}{$t$}\!\left/ \!\raisebox{-1ex}{$t_R$}\right.\right)}^{n-1}}{\left(n-1\right)!}{F}_n\left(x,z,D,H\right)$$

(4)

The surface velocity of the viscoelastic model is related to the elastic layer thickness (H), the earthquake recurrence period T, the relaxation time τ_M, and other parameters [109]. A prior constraint of the unknown parameters is employed according to

Table A1. A Priori Bounds in inversion for ATF.

Parameter	Minimum	Maximum	Priori Value
Fault slip rate (mm/yr)	−20	20	10
Locking depth (km)	0	35	18
Velocity offset (mm/yr)	0	15	3
Recurrence interval (years)	800	1500	1000
Time since last earthquake (years)	100	100	100
Elastic thickness of FNE, H1 (km)	10	65	35
Elastic thickness of FSW, H2 (km)	10	65	20
Viscosity η	1016	1022	1019
Viscoelastic lower crust η2	1016	1022	1019
Rigidity ratio, μ1/μ2	0	5	2

.

For a ground point with z = 0, the $F_n$ is simplified as [107]:

$$F_n\left(x,z=0,D,H\right)=ta{n}^{-1}\left(\frac{x}{2nH-D}\right)-ta{n}^{-1}\left(\frac{x}{2nH+D}\right)$$

(5)

In the case of model B, the best-fitting model parameters are shown in Figure 7b and Figure 8b. For the Nierkule segment (MN in Figure 4), the inferred long-term slip rate is 7.3 ± 1.1 mm/yr, the locking depth is 15.8 ± 4.3 km, and the elastic layer thickness is 26.1 ± 2.2 km. For the Kurukuole section (PQ in Figure 4), the optimal values are estimated to be 7.5 ± 0.3 mm/yr slip rate, 14.8 ± 4.9 km locking depth, and 19.8 ± 4.3 km elastic layer thickness. The slip rate obtained by this model is slightly larger than that obtained by the buried screw dislocation model in elastic half-space, but it is still smaller than the long-term geological rate.

Model C is still a viscoelastic model adopted in [107]. The difference is that different effective elastic layer thicknesses on both sides of the fault are considered. To reduce the number of parameters in the inversion, the locking depth is fixed as the results from Model B and no longer participates in inversion because of the trade-off between fault slip rate and locking depth. We assume that the Tarim Basin and the TP have different effective elastic thickness and shear modulus, respectively, and other prior conditions are the same as Model B to inverse the difference in elastic layer thickness and shear modulus ratio on both sides of the fault.

The inversion results are shown in Figure 7c and Figure 8c. The long-term slip rate of the Nierkule segment (MN profile inversion) estimated by the model is 9.8 ± 1.1 mm/yr, the viscosity coefficient of the viscoelastic layer in the lower crust is 3.2 ± 3.3 × 10²⁰ Pa s, and the thickness of the elastic layer on the side of the Tarim Basin is 36.3 ± 6.1 km. The thickness of the elastic layer on the side of the Tibetan Plateau is 18.3 ± 2.3 km, and the ratio of the shear modulus on both sides is 1.8 ± 1.2. For the Kurukuole segment (PQ profile), the inferred slip rate is 8.6 ± 1.1 mm/yr, and the viscosity coefficient in the lower crust is 1.7 ± 2 × 10²⁰ Pa s. The thickness of the elastic layer on the side of the Tarim Basin is 34.8 ± 4.2 km, and that on the other side of TP is 18.8 ± 2.9 km. The ratio of shear modulus on both sides is 1.5 ± 1.1, and the inversion results of fault slip rate by this model are close to the long-term geological rate. The best-fitting model is shown in Figure 9.

In addition, two other cases are also tested. One is that the thickness of the elastic layer on both sides of the fault is the same, but the shear modulus is different (that is, the ratio parameter of shear modulus on both sides of the fault is added based on Model B). The test results showed that the long-term slip rate increased to 9.5 ± 0.5 mm/yr and 9.3 ± 0.5 mm/yr, respectively. The ratio of the shear modulus of the upper crust on the side of the Tarim Basin to that on the other side of TP is 0.5 ± 0.2 (see Figure A2), which means that the rigidity of the upper crust on the side of Tarim Basin is weaker than that on the other side of TP, which is not consistent with general understanding. In the other case, the shear modulus on both sides of the fault is the same, but the elastic layer thickness is different (the shear modulus ratio parameter is removed based on Model C). The test results are shown in Figure A3. The inferred long-term slip rate estimated is close to model C, which is 9.3 ± 0.2 mm/yr and 9.4 ± 0.2 mm/yr, respectively. However, the obtained elastic layer thickness on the side of Tarim Basin is large, the inversion result for the MN profile is 41.5 ± 2.5 km, and the inversion result for the PQ profile is 50.6 ± 3.4 km. The elastic layer thickness on the side of TP does not change much with Model C. That is to say, the big difference in elastic layer thickness on both sides of the fault can also cause such asymmetry of deformation field in the Altyn fault zone. The study of Te based on topographic and gravity data shows that the Tarim Basin has a moderate elastic layer thickness, with Te values ranging from 20 to 60 km. The closer it is to the TP, the lower the value is, and the higher value, the near center of the basin is [110,111,112]. Intermediate Te values appear in the Qaidam Basin, filled with thick sediments. The Altyn and East Kunlun regions obviously have low elastic layer thickness, with Te values ranging from 10 to 30 km. Generally, large deformation also occurs in these weak zones, which are also the main areas absorbing the shortening of the India-Euroasia continental collision [112]. It can be considered that the addition of the shear modulus ratio parameter does not affect model inversion. As a characterization of lithospheric strength, the spatial variation of effective elastic thickness is of great significance to lithospheric deformation. In addition, because elastic thickness mainly depends on crustal thickness, temperature, composition, etc., it can be used to reflect the lateral changes of deep lithosphere structure, and the elastic layer thickness on both sides of the fault has a great influence on the inversion results. Results based on the four combined test results of models B and C and the previous research results on the thickness of the elastic layer around the TP, we prefer the results of Model C. The differences in elastic layer thickness and shear modulus on both sides of the fault jointly cause the asymmetry in deformation across the ATF.

4.1.3. Three-Layer Model

The large-scale tectonic evolution of the Qinghai Tibet Plateau and its marginal areas show that the weak ductile lower crust plays an important role. Therefore, Model D, following the practice of DeVries and Meade [113], establishes a simple three-layer viscoelastic model on the side of the TP to restrict the viscosity and thickness of the lower crust of the ATF, a low viscosity viscoelastic channel (H1) with a thickness of 20 km is embedded on the right side between the elastic block and the bottom substrate. The two-layer elastic, viscoelastic model (as shown in Figure 6d) is still used on the side of the Tarim Basin. The inversion results are shown in Figure 7d and Figure 8d. The inferred long-term slip rate of the Nierkule segment (MN profile) is 7.2 ± 0.2 mm/yr, and the viscosity coefficient of the lower crust is 2.9 ± 2.5 × 10²⁰ Pa s, the viscosity coefficient of the embedded layer is 4.7 ± 4.1 × 10¹⁹ Pa s, the thickness of the elastic layer on the side of the Tarim Basin is 47.6 ± 4.5 km, and on TP side is 15.8 ± 6.5 km. For the Kulukuole section (PQ profile), the inferred long-term slip rate is 6.9 ± 0.2 mm/yr, and the viscosity coefficient of the lower crust is 0.9 ± 0.6 × 10²⁰ Pa s, the viscosity coefficient of the embedded layer is 3.0 ± 0.6 × 10¹⁹ Pa s, the thickness of the elastic layer on the Tarim Basin side is 37.4 ± 4 km, and the thickness of the elastic layer on the TP side is 9.1 ± 0.7 km [114]. The geomagnetic data revealed the existence of two huge middle and lower crustal low-resistance anomalous zones, which are considered as two weak material flows in the middle and lower crust, where thermal fluids are abundant [115]. The purpose of Model D is to test whether there is a low viscosity viscoelastic channel on one side of the TP. Compared with the Altyn collision orogenic belt where the Model D section is located, the model may be more suitable for the Lhasa block and Qiangtang terrane suture belt.

4.1.4. Summary

In short, Model A is a half-space model with different rigidity on both sides of the fault, highlighting the asymmetric pattern is only due to a lateral variation in rigidity. Models B and C are viscoelastic-coupling models, where both the elastic thickness and shear modulus are free parameters, highlighting the difference between elastic layer thickness and shear modulus on both sides of the fault might contribute to the asymmetric strain pattern across the western ATF. The Model D is a three-layer model in which a viscoelastic mid-crustal layer is sandwiched between an elastic upper crust and an elastic upper mantle, aiming to test whether the middle to lower crust may be weak.

Integrating the above comparisons and results, we conclude that the difference in elastic layer thickness and shear modulus might lead to the asymmetric strain across the ATF. Nevertheless, our current study could not distinguish which factor plays a dominant role. Further studies are needed in the future.

In conclusion, there are differences in the physical properties of the media on the north and south sides of the ATF, and the above results are verified by inversion. This difference may be related to the collision orogeny between the Tarim block and the Altyn block. In addition, there are some shortcomings in this paper. On the north side of the ATF, there is a parallel Cherchen fault, which should be classified into the Altyn fault zone System from the perspective of its distribution characteristics and kinematic characteristics. On the south side, there are also several large faults whose deformation characteristics are shown in the InSAR deformation field. The large-scale surface deformation signal observed by satellites should be the interaction of these faults. Here, it is only simplified to be controlled by a single fault, and the research of multi fault combination model can be strengthened in the future.

4.2. Interseismic Slip Rate

As a boundary strike-slip fault, the fault slip rate of the ATF is useful to constrain the current continental dynamic process [116], as well as to further explore the deformation mechanism of the TP and quantitatively evaluate seismic risk. Therefore, a great deal of research work has been carried out to obtain the slip rate of each segment of the ATF by using Geodesy (GPS, InSAR) and geological methods.

The geological method is mainly based on accumulated displacements of the offset landforms and the corresponding dating data to estimate the average slip rate of fault movement. Xu et al. [117] and Meriaux et al. [6,118] obtained a slip rate of 17–20 mm/yr through a detailed study of each segment of the ATF. Zhang et al. [119] and Cowgill et al. [120] resummarized and analyzed previous research data and believed that the ATF slip rate should be ~10 mm/yr. Cowgill et al. [120] assumed that the abandonment age for the upper terrace was the offset time of the alluvial terrace scarp, and the calculated geological slip rate was 9.4 ± 2.3 mm/yr. Zhang et al. [119] also recalculated the fault slip rate from the abandoned age of upper and lower terraces and the activity of fault scarps and yielded a result consistent with the slip rate observed by GPS and InSAR, which could solve the mismatch between the geological slip rate and the slip rate obtained by GPS. Later, Gold et al. [121] and Cowgill et al. [122] obtained geological slip rates of different sites along with the ATF, which were basically at low slip rates. In the past 20 years, with the advance in dating of the geomorphological markers, it seems that both the slip rate obtained by geodesy [20,22,123,124,125] and the slip rate obtained by geology in recent years [123,125,126] have been well consistent. That is, along most of the ATF, there is a uniform slip rate of ~10 mm/yr. Our results show that after using the layered model, considering the vertical variation of viscosity coefficient, the long-term slip rate is 9.8 ± 1.1 mm/yr, and the estimated slip rate of ATF is consistent with the geological data. If the viscous effect of the lower crust and upper mantle is ignored, it will lead to a misestimation of the slip rate.

5. Conclusions

The fault slip rate and locking status of the boundary faults are important for understanding regional tectonic deformation and seismogenic potential. The rapid development of spatial geodesy techniques (GPS and InSAR) provides technical support for obtaining high-resolution and high-precision deformation maps, and studying fault slip rate and locking depth during the earthquake cycle. Our key conclusions are summarized as follows:

• Different models are used to invert the fault slip rate and locking depth of the western segment of ATF. In the viscoelastic model, the media on both sides of the fault are divided into the elastic upper crust, viscoelastic lower crust, and upper mantle, and stratified lateral inhomogeneity parameters are considered. From the west to east section, the inferred sinistral strike-slip rate of the ATF is 9.8 ± 1.1 mm/yr and 8.6 ± 1.1 mm/yr, respectively, and the locking depth is 15.8 ± 4.3 km and 14.8 ± 4.9 km respectively. Based on the elastic model fitting, the left-lateral strike-slip rates of the ATF are 7.1 ± 0.1 mm/yr and 6.1 ± 0 mm/yr, respectively, and the locking depths are 18.1 ± 0.6 km and 11.7 ± 0.4 km, respectively. The results show that ignoring the viscoelastic effect will significantly affect the estimation of fault slip rate and locking degree;
• The decomposition of the ascending and descending LOS velocities into fault parallel and vertical components in this paper shows that although the ATF is mainly strike-slip, it is accompanied by a small amount of thrust, and the thrust component gradually increases from west to east, reaching 1 mm/yr in the Altyn Mountains. If the vertical deformation of InSAR data is ignored, it may bias the interpretation of crustal deformation. Constrained by the fault parallel velocity field from west to east, the sinistral strike-slip rates of the ATF obtained in this paper are 9.8 ± 1.1 mm/yr and 8.6 ± 1.1 mm/yr, respectively, in line with the characteristics that the ATF is in low slip rates and decreases gradually from west to east which is obtained by GPS inversion;
• The velocity profile across the ATF shows that there is asymmetry across the fault. The results integrate the combined test results of different parameters of Models B and C and the previous research results on the thickness of the elastic layer around the TP; we believe that the difference between elastic layer thickness and shear modulus on both sides of the fault jointly causes the asymmetry of interseismic velocity.