Rheology of the Northern Tibetan Plateau Lithosphere Inferred from the Post-Seismic Deformation Resulting from the 2001 Mw 7.8 Kokoxili Earthquake

Abstract

Time-series observations from Sentinel-1 A/B InSAR spanning November 2014 to July 2021 were used to study the late post-seismic deformation velocity field arising from the Kokoxili earthquake. The deformation velocity caused by the interseismic slip along the major active faults in Tibet was first simulated. Comparing the simulated deformation velocity with the observed one, the maximum ratio of the simulated deformation velocity to the observed one was found to be 42%, indicating continuity in the viscoelastic relaxation caused by the 2001 Kokoxili earthquake. Subsequently, the rheological structure of the Kokoxili region was explored using a mixed model comprising the viscoelastic relaxation mechanism and the buried elastic dislocation model. The best estimated viscosities for the lower crust and upper mantle were ${\mathsf{\eta }}_{\mathrm{lc}}=1_{-0.44}^{+0.78}\times {10}^{19} \mathrm{Pas}$ and ${\mathsf{\eta }}_{\mathrm{um}}=1_{+0}^{+0.78}\times {10}^{20} \mathrm{Pas}$, respectively. The results obtained in this study were compared with those of previous studies that used the early post-seismic displacement ranging from 0 to 6.5 years following the earthquake. The obtained value was largely the same as the previously estimated steady-state viscosity, which means that the viscosities of the viscoelastic layer beneath the Kokoxili regions have almost reached their stable state. Furthermore, the effective lower crustal viscosity of the Kokoxili region exhibited a logarithmic trend with time.

1. Introduction

The Mw 7.8 Kokoxili earthquake of November 2001 occurred along the Kunlun fault and was located on the Songpan–Ganzi terrane [1,2] of the Northern Tibetan Plateau. It was the largest to have occurred along the Kunlun fault over the past 100 years [3] and had a significant influence on the stress field of Northern Tibet. The 2001 Kokoxili earthquake provides an unprecedented opportunity to investigate the rheological properties of the lithosphere in this part of Tibet using observations of post-seismic deformation.

The main mechanisms proposed for assessing the post-seismic deformation include (i) afterslip within or near the coseismic rupture area [4,5,6,7]; (ii) viscoelastic relaxation of the ductile layers, where coseismically induced stress variations in the lower crust and upper mantle cannot be sustained and are released through viscoelastic flow [8,9,10]; (iii) poroelastic rebound, where the fluids in the upper crust flow due to the coseismic pressure changes [11,12,13]. These three mechanisms have different spatial and temporal behaviors. Poroelastic rebound and afterslip mainly occur in the near field and can last for a few months to years. Viscoelastic relaxation mainly occurs in the far-field and can last for decades. Identifying the mechanics responsible for the post-seismic deformation can help provide insight into future seismic hazards because the released post-seismic stress may push adjacent segments closer to failure [14,15,16].

The viscoelastic relaxation process provides a convincing tool to investigate the strength of the continental lithosphere which is causing considerable controversy [17]. By exploring the viscoelastic relaxation following large earthquakes, bounds can be placed on the strength of the lower crust and upper mantle. Previous studies have shown that viscoelastic relaxation can occur in the lower crust in the case of some earthquakes, in the upper mantle in other cases, and sometimes in both. Lower crustal relaxation has been observed in the 1980 Mw 6.9 Italy Irpinia earthquake [18], 1989 Mw 7.1 Loma Prieta earthquake [19], 1990 Mw 6.4 Gonghe earthquake [20], 1994 Mw 6.7 Northridge earthquake [21], 1997 Mw 5.6–6 Italy Umbria Marche earthquake sequence [22], 1997 Mw 7.6 Manyi earthquake [23], 1999 Mw 7.6 Chi-chi earthquake [24], 2001 Mw 7.8 Kokoxili earthquake [2], 2001 Mw 7.6 Bhuj earthquake [25], 2003 Mw 6.6 Bam earthquake [26], 2004 Mw 6.0 Parkfield earthquake [27], and 2008 Mw 7.8 Wenchuan earthquake [28], with observed data spanning 3 months to 7 years following the earthquakes. In these regions, the viscosity of the lower crust is lower than that of the upper mantle, indicating that the lower crust is weaker than the upper mantle. Earthquakes with viscoelastic relaxation processes occurring in the upper mantle are mainly located in the western US, such as in the 1915–1954 earthquake sequences occurring along the Central Nevada Seismic Belt [29], 1959 Mw 7.3 Hebgen Lake earthquake [30], 1992 Mw 7.3 Landers earthquake [31], 1999 Mw 7.1 Hector Mine earthquake [32], and 2002 Mw 7.9 Denali earthquake [33] with the post-seismic time period ranging from the first 0–9 months to 0–27 years. Upper mantle relaxation has also been found in the 2000 Mw 6.5 Iceland earthquakes [34] and 2010 Mw 7.2 El Mayor Cucapah earthquakes [35] using the first 0–5 years and 0–1.5 years post-seismic data, respectively. In these regions, the viscosity of the upper mantle is lower than that of the lower crust, which means that the upper mantle is weaker than the lower crust. Only a few studies have found the relaxation phenomenon in both the lower crust and upper mantle; representative earthquakes include the 1999 Mw 7.5 Izmit Turkey earthquake [36] and the seven Mongolia earthquakes [37] that hit during the 1900s with magnitudes ranging from Mw 6.8 to 8.4; the studied post-seismic time periods were the first 0–6 years and 1997–2004, respectively. Furthermore, most studies focused on early and short-term viscoelastic relaxation following the earthquakes, with few on late and long-term viscoelastic relaxation. Long-term viscoelastic relaxation has important research value. For example, Gualandi [38] pointed out that the long-term relaxation following the 2010 Mw 7.2 El Mayor-Cucapah earthquake affected the interseismic velocities and regional seismicity rates for more than a decade. Fujiwara [39] studied the temporal function of post-seismic deformation following the 2011 Tohoku–Oki earthquake using data from more than 10 years and found that the long-time function can help understand the physical phenomena. Therefore, it is meaningful to conduct a long-term relaxation study of earthquakes.

In the case of the Kokoxili earthquake, although there have been many post-seismic deformation studies [2,3,23,28,40,41,42], early observations of the post-seismic deformation were used to estimate the viscosities, and the longest period considered was the first 6.5 years after the event [3]. Whether the viscoelastic relaxation process continues to this day remains unclear. Yanchuan [43] believed that post-seismic deformation may have continued up to 2014 based on the time trend in the post-seismic deformation. However, no quantitative analysis was performed. Furthermore, the viscosity following an early post-seismic deformation is often lower than that following a late post-seismic deformation [28]. Thus, it would be meaningful to investigate whether the post-seismic deformation continues to exist and explore the rheological structure using observations of the post-seismic deformation in this period.

In this study, observations from the Sentinel-1 mission spanning November 2014 to July 2021 (late post-seismic deformation, 13–20 years after the Kokoxili earthquake) were used to investigate the post-seismic deformation velocity following the Kokoxili earthquake. The observed deformation velocity was first compared with the simulated one caused by the interseismic slip along the major active faults in Tibet, and the effective viscosities in this late post-seismic period were estimated using a mixed model comprising the viscoelastic relaxation mechanism and the buried elastic dislocation model. Finally, the best estimate viscosities were compared with the results of previous studies that used early post-seismic deformation, and whether the viscosities recovered to a relatively stable state was investigated.

2. Geologic Setting

The Tibetan Plateau was developed following the Indian continental subduction and collision with the Eurasian plate some 50 Myr ago, and the convergence rate is approximately 3 cm/yr [44]. The collision caused the thickening of the Tibet Plateau crust, and the crust thickened from approximately 55 km in the northern plateau to approximately 80 km in the southern plateau [2]. Due to the collision, large strike-slip faults were formed in Tibet and thrust faults at the boundary of Tibet. The major active strike-slip faults are the Altyn Tagh fault, Manyi fault, Kunlun fault, Haiyuan fault, Yushu–Ganzi fault, and Xianshuihe fault (Figure 1), and the largest thrust fault is the Himalaya Frontal Thrust (HFT). The Tibetan Plateau has been undergoing intense tectonic deformation; however, the formation of the Tibetan Plateau and its current deformation status remain outstanding questions and are under debate [45]. The key to answering these questions is to understand the rheological nature of the lithosphere beneath the plateau. Two classes of models have been proposed to explain the growth of the Tibetan Plateau: (i) crustal rigid blocks [46,47] causing lateral extrusion along the strike-slip faults, and (ii) channelized lower crustal flow [48,49]. For the second model, it remains unclear how much of the lithosphere is weak enough to flow and at what depth [2] for the different parts of the Tibetan Plateau.

Tibet is frequently hit by strong earthquakes (Figure 1), which are located at the faults in or around the plateau. North Tibet has witnessed several strong earthquakes, such as the 1997 Manyi earthquake, the 2001 Kokoxili earthquake, and the 2021 Maduo earthquake. North Tibet can be mainly divided into three parts: Qiangtang terrane, Songpan–Ganzi terrane, and Qaidam Basin [28]. Previous studies on the rheological structure of North Tibet often divided the lithosphere into the elastic upper crust, viscoelastic lower crust, and viscoelastic upper mantle. The depth of the elastic upper crust is approximately 15–20 km where most earthquakes have occurred [23]. The depth of Moho is approximately 70 km based on the observation of seismic wave data and gravity data [50,51,52]. The strength of the lower crust is equal to or lower than that of the upper mantle based on the high electrical conductivity [53], lower shear save speed [54], high Poisson ratio [55], seismic anisotropy variations along the depth [56], and post-seismic deformation analyses [2,3,23,28,40,41,42]. Furthermore, the Kunlun and Altyn Tagh faults are two major strike-slip faults in northern Tibet. They accommodate the eastward extrusion of Tibet in response to the Indo-Asian collision. Among them, the 1500 km-long east–west striking Kunlun fault (Figure 1) is a left-lateral strike-slip fault. It is located in the Kunlun Mountain range and accommodates the eastward motion of the Qiangtang and Songpan–Ganzi terranes relative to the Qaidam Basin [22]. The slip rates along most of the Kunlun fault has been estimated to be approximately 10–12 mm/yr [2,22] based on geological and GPS data. This relatively high slip rate means that the Kunlun fault plays a key role in the extrusion model for northern Tibet. The Kunlun fault can be broadly divided into five main segments, with the slip rate decreasing from west to east segments [43]). Five large left-lateral strike-slip earthquakes (M ≥ 7) have occurred along the Kunlun fault over the last century [2]: the 1937 Ms 7.5 Huashi Canyon earthquake, the 1963 Ms 7.1 Dulan earthquake, the 1973 Mw 7.4 Manyi earthquake, the 1997 Mw 7.6 Manyi earthquake, and the 2001 Mw 7.8 Kokoxili earthquake. The Huashi Canyon and Dulan earthquakes ruptured the central segment of the Kunlun fault, with the epicenter moving increasingly west [41]. The 1973 and 1997 Manyi earthquakes both ruptured the western part of the Manyi fault, which is one of the splays split from the Kunlun fault [57]. The 2001 Kokoxili earthquake hit along the westernmost region of the Kunlun fault.

The epicenter of the 2001 Kokoxili earthquake determined by the US Geological Survey (USGS) was 35.95N, 90.54E. The magnitude of this left-lateral slip earthquake was Mw 7.8. Geodetic and seismic observations have offered valuable constraints on coseismic slip distribution models describing this event [58,59,60]. These studies have concluded that this earthquake occurred along a vertical fault and that coseismic slip occurred between depths of 0 and 16.5 km. In this study, the coseismic slip model proposed by Wang [58] was employed for the viscoelastic relaxation simulations.

3. InSAR Data and Processing Methodology

The Interferometric Wide swath (IW) mode of Sentinel-1 (

Table 1. Sentinel-1 data list used in this work.

Path	Time Period	Number of Images
AT12	2014/11–2021/07	165
AT114	2014/11–2021/05	159
AT41	2014/11–2021/07	163
AT143	2014/11–2021/05	155
AT70	2014/11–2021/07	158
AT172	2014/11–2021/07	161
DT19	2014/11–2021/07	161
DT121	2014/11–2021/07	158
DT48	2014/11–2021/03	151
DT150	2014/11–2021/05	154
DT77	2014/11–2021/07	143

) was used in this study. A total of 961 ascending Sentinel-1 A/B images from Path 12, 114, 41, 143, 70 and 172, and 767 descending Sentinel-1 A/B images from Path 19, 121, 48, 150, and 77 (from November 2014 to July 2021 for both Paths) were used. Note, that there are data gaps of 0.25 years and 0.3 years for ascending Paths 114 and 143 and descending Paths 48 and 150, respectively. Given that they represent only 3% and 4% of the 6.8-year study period, this data scarcity is not a problem.

Interferograms with seven looks in the azimuth direction and 21 looks in the range direction were constructed using the JPL/Caltech ISCE stack processor [61,62]. For each epoch, we generated interferograms using itself and the next five epochs following it. The 30 m Shuttle Radar Topography Mission (SRTM) digital elevation model [63] and the Precise Orbit data were used to simulate and remove the topographic phases and flatten earth phases from each interferogram. By finding the offsets between the master SLC (Single Look Complex) and SLCs using the DEM and orbit vectors, multi-looked and filtered interferograms were co-registered to a single master SAR image. Finally, the statistical-cost network-flow algorithm (SNAPHU) [64] was used to unwrap the phases of the co-registered interferograms. Based on these steps, we obtained a stack of phase-unwrapped interferograms co-registered to a common SAR acquisition, corrected for earth curvature and topography.

The open-source Miami InSAR time-series software in Python (MintPy) is used for time-series processing [65]. MintPy uses distributed scatterers and is an improved small baseline subsets (SBAS) algorithm, which uses a fully connected network of interferograms and performs phase corrections in the time-series domain, in contrast to the conventional interferogram domain. Based on the routine processing workflow of MintPy, the raw phase time-series was first inverted from the interferogram network (known as phase triangulation). Subsequently, the solid earthquake tide phase, topographic residuals, and tropospheric delay using the global atmospheric models (ERA5 from the ECMWF) were calculated and removed from the raw phase time-series. Finally, the average line-of-sight (LOS) velocity was estimated on a pixel-by-pixel basis using the noise-reduced displacement time series.

For each pixel, the deformation is a relative measurement with respect to a reference pixel on the same track. Thus, when the deformation velocity fields of adjacent tracks are concatenated, a constant offset, which is the median [66] of the differences in the deformation velocity fields in the overlapping areas, is first estimated. Subsequently, the average deformation velocity of the two tracks for the overlapping areas is used.

Finally, the ascending and descending data are combined to vertical and east-west deformation velocity components, assuming that the ascending images are acquired on the same day as the descending images, although they are acquired four days later. Ascending data with linear ramps removed are presented throughout the paper. Liang [67] found that the diurnal variations of the ionosphere make Sentinel-1 ascending data have much stronger ionospheric effects (Figure S1).

4. Modelling Approach

4.1. Interseismic Deformation Velocity Field Simulation

In this study, the interseismic deformation velocity field caused by the aseismic slip along the major active faults in Tibet was simulated. The faults are Altyn Tagh fault, Manyi fault, Kunlun fault, Haiyuan fault, Yushu–Ganzi fault and Xianshuihe fault (Figure 1). A buried elastic dislocation model in an elastic half space [68] was used to calculate the aseismic surface velocity relative to the deep fault slip movement [69,70]. For all faults except for the Kunlun fault, one fault dislocation from the locking depth to the depth of 1000 km (representing the infinite elastic half space) was used for the simulation (

Table 2. Locking depth and average slip rate of faults.

Fault	Length [km]	Locking Depth [km]	Average Slip Rate [mm/yr]
Alytn Tagh [71]	2000	20	7
Manyi [72]	200	22	5
Haiyuan [73]	670	15	4.5
Yushu-Ganzi [70]	668	21	15
Xianshuihe [74]	379	20	9.2

). The Kunlun fault was into five segments with different slip rates (

Table 3. Locking depth and slip rate of Kunlun fault segments [43,75].

Fault Segment	Length [km]	Locking Depth [km]	Slip Rate [mm/yr]
Kusai Lake	283	20	15
Xidatan–Dongdatan	135	20	17.4
Alake Lake	164	20	10.7
Tuosuo Lake	224	20	6.2
Maqin Maqu	359	20	6.1

); the 2001 Kokoxili earthquake ruptured the entire Kusai Lake segment and the western part of the XidatanvDongdatan segment.

4.2. Setting of Viscoelastic Relaxation Model Parameters

In this study, the open-source program RELAX [76,77] was used to calculate the viscoelastic relaxation velocity and the coseismic parameters are obtained from the slip model of Wang [59] for the Kokoxili earthquake.

When simulating viscoelastic relaxation deformation velocity, a Maxwell rheology, which has two parameters (steady-state viscosity and steady-state shear modulus), was used for the two viscoelastic layers whose parameters need to be estimated, namely the lower crust and the upper mantle. To model the viscoelastic relaxation with a Maxwell rheology, a three-layered rheology model was used to investigate a range of viscosities. The elastic parameters of the models were kept constant throughout. The shear modulus was set to 30 GPa and the Poisson ratio was 0.25. The upper layer is the elastic upper crust whose thickness is fixed at 16.5 km, which is consistent with the depth above which most coseismic slip occurred [58,59,78]. The second layer is the viscoelastic lower crust with a viscosity (${\mathsf{\eta }}_{\mathrm{lc}})$ ranging from $1\times {10}^{17} \mathrm{Pas}$ to $1\times {10}^{21} \mathrm{Pas}$. From CRUST1.0 [79] and the observation of teleseismic shear-coupled P waves [55], it can be inferred that the Moho depth of the Songpan–Ganzi terrane is approximately 55km. However, the observation of broadband teleseismic body waves [50], the gravity data [51] and the wide-angle reflection and refraction profile [52] collected across the Northeast Tibet Plateau show the crustal thickness of Songpan–Ganzi terrane and Kunlun Mountain to be approximately 70 km. Moreover, Ryder and Yangmao both used a crustal thickness of 70 km [2,41]. Therefore, the base of the crust is fixed at 70 km in this study. The third layer is a viscoelastic half-space corresponding to the upper mantle and whose viscosity (${\mathsf{\eta }}_{\mathrm{um}}$) ranges from $1\times {10}^{17} \mathrm{Pas}$ to $1\times {10}^{21} \mathrm{Pas}$. The step size of the considered viscosities is ${\log }_{10}\mathsf{\eta }=0.25$.

4.3. Inversion Approach

A mixed model was used to estimate the best-fitting viscosity. The mixed model combines the viscoelastic relaxation mechanism with the buried elastic dislocation model. In the buried elastic dislocation model, the surface deformation is modeled by the deep slip of a simple dislocation locked from the locking depth to the surface and slipping below to an infinite depth in an elastic half space [68]. Since a precise slip rate distribution along the Kunlun fault is beyond the scope of this paper, the Kunlun fault is simplified as one finite dislocation with a locking depth of 20 km and is allowed to accommodate the slip rate of other major active faults. Besides the depth, the strike and dip are fixed at 100° and 90°, respectively, and the length, strike slip and dip slip were inverted. For the viscoelastic relaxation model, we inverted the viscosities of the lower crust and upper mantle.

The uniform downsampling method was used to sample from the observation (see Supplementary Materials Figures S2 and S3 for downsampling results). The quadtree downsampling method was not used because its results cannot capture the detailed deformation north of the fault (Figure S4). The misfit function is calculated as follows:

$${\chi }^2={\left({\mathit{d}}_{\mathit{o}\mathit{b}\mathit{s}}-{\mathit{d}}_{\mathit{s}\mathit{i}\mathit{m}}\right)}^T{\mathit{C}}^{-1}\left({\mathit{d}}_{\mathit{o}\mathit{b}\mathit{s}}-{\mathit{d}}_{\mathit{s}\mathit{i}\mathit{m}}\right)$$

where ${\mathit{d}}_{\mathit{s}\mathit{i}\mathit{m}}$ is the simulated deformation velocity, ${\mathit{d}}_{\mathit{o}\mathit{b}\mathit{s}}$ is the observed deformation velocity, and $\mathit{C}$ is the covariance matrix, calculated following Bagnardi and Hooper [80]. ${\mathit{d}}_{\mathit{o}\mathit{b}\mathit{s}}$ contains 920, 944, 1094 observations and 926, 1012, 948 observations for ascending track 114 & 41, 143, 70 & 172 and descending track 121 & 48, 150 and 77, respectively. The data were sampled from each of the AT 143, DT 150 and DT 77 swaths without concatenation i.e., without averaging the overlap areas. Concatenation was performed for AT114 with AT41, AT70 with AT172, and DT 121 with DT48 and treated as one swath. This is because the studied deformation area (red rectangle in Figure 1) only covers a small part of AT114, AT172 and DT48. For each of the swaths, we invert for linear ramps (constant offset, azimuth ramp and range ramp) [80].

Figure 2 shows the flowchart of the inversion approach.

First, the RELAX program [76,77] was used to calculate the predicted east–west, north–south, and vertical direction velocities due to the viscoelastic relaxation for the given viscosities. The predicted velocity in the three directions was projected onto the radar LOS of the ascending and descending viewing geometries.

Second, these predicted LOS relaxation velocities are subtracted from the InSAR observed LOS velocities to obtain the observation residual deformation velocities.

Third, a Bayesian approach with the GBIS software was used to find the additional deep slip rate along the Kunlun fault.

Finally, the misfit was obtained for the given viscosities in step 1 and the best-fitting buried elastic dislocation. Finally, the best estimate viscosities correspond to the minimum misfit.

5. Results

5.1. Observed Post-Seismic Deformation Velocity

Figure 3 shows the concatenated LOS deformation velocity over the period November 2014 to July 2021. The ascending data show large lobes of LOS increase and LOS decrease velocity (up to 0.7 cm/yr and 1.0 cm/yr, respectively) south and north of surface rupture of the Kokoxili earthquake. The descending data show a narrow rectangle LOS decrease (up to 0.8 cm/yr) north of surface rupture and a discrete lobe of LOS increase (up to 0.6 cm/yr) south of surface rupture.

There are notable seams bordering the overlap areas in the concatenated data. These seams can be attributed to two factors. One is that the median of the velocity differences is used as constant offset, which contains errors because of the differences of the velocity difference from the median value; the other is that the errors of the two concatenating tracks are different.

5.2. Simulated Interseismic Deformation Velocity

Considering that the deformation velocity of the late post-seismic period is used in this study, the magnitude of the interseismic deformation signal should be discussed first. We mainly compared the simulated and observed east–west deformation velocity because the major active faults in Tibet are strike-slip faults. Figure 4b shows the simulated east-west deformation velocity, in comparison with the observed east-west deformation velocity (Figure 4a). The deformation spatial patterns of the observation and simulation data in the Kokoxili region (black rectangle in Figure 4) are similar. There is a large lobe with significant deformation south of the Kokoxili fault, and the deformation region north of the Kokoxili fault is relatively narrow.

The maximum simulated east deformation velocity is 0.75 cm/yr, accounting for approximately 42% of the maximum observed east deformation velocity of 1.8 cm/yr, and the maximum simulated west deformation velocity is 0.2 cm/yr, accounting for approximately 13% of the maximum observed west deformation velocity of 1.5 cm/yr. Therefore, it is reasonable to speculate that there is still viscoelastic relaxation following the Kokoxili earthquake from November 2014 to July 2021.

5.3. Inversion Results of the Mixed Model

The regional layer structure model, the viscosity range, and step size considered for the lower crust and upper mantle in the viscoelastic relaxation model are described in Section 4.2. We fixed the width, depth, strike and dip of the buried elastic dislocation model and inverted for the length, strike slip and dip sip.

Our results (Figure 5b) yield the best-fitting parameters ${\mathsf{\eta }}_{\mathrm{lc}}=1_{-0.44}^{+0.78}\times {10}^{19} \mathrm{Pas}$ and ${\mathsf{\eta }}_{\mathrm{um}}=1_{+0}^{+0.78}\times {10}^{20 }\mathrm{Pas}$. The minimum misfit is obtained at ${\mathsf{\eta }}_{\mathrm{lc}}=1.78\times {10}^{19} \mathrm{Pas}$ and ${\mathsf{\eta }}_{\mathrm{um}}=1\times {10}^{20} \mathrm{Pas}$. The best estimates of length, strike slip rate and dip slip rate are 215 km, 17 mm/yr and 3 mm/yr, respectively (

Table 4. Range, optimal and uncertainties of InSAR-inversion fault slip parameters for the buried elastic dislocation model.

	Length (km)	Width (km)	Depth (km)	Dip (°)	Strike (°)	Strike Slip (m/yr)	Dip Slip (m/yr)	Lat a	Lon a
Lower	20	980	20	90	100	−0.1	−0.1	—	—
Upper	1200	Same as above				0.1	0.1	—	—
Optimal	215.17	—	—	—	—	0.017	−0.003	35.85	92.83
2.5%	102.76	—	—	—	—	0.006	−0.02	—	—
97.5%	305.26	—	—	—	—	0.023	0.01	—	—

^a longitude and latitude of the middle point of fault upper edge.

, Figure 6).

A comparison between the observed and simulated LOS deformation velocities of each track (see Supplementary Materials Figures S5 and S6) shows that the model can explain the descending and ascending data well. For a more intuitive presentation of the data-model fitting results, the ascending and descending tracks for the data, model and residual are concatenated. The concatenated results (Figure 6) once again show that the model can explain the observation well, particularly for the large deformation lobe south of the surface rupture.

However, there are two small deformation areas, the northwest corner of the fault in ascending data and the southeast corner of the fault in descending data, which cannot be explained well by the best estimated model. The possible reasons are 1. A single buried, uniform dislocation was used for modeling, and it was assumed that the entire seismogenic zone [0–20 km] was fully coupled; 2. The Kunlun fault was used to accommodate the observed interseismic velocity generated by the major active faults in Tibet. This simplified fault model may not only lead to a relatively overestimated slip rate but may also fail to explain the LOS increase signal at the southeast corner of the fault in the descending data. This signal could be caused by the Ganzi–Yush fault.

Although there are two small deformation areas that cannot be explained well by the best estimated model, the best-fitting model can explain most of the deformation on both sides of the Kokoxili fault well. Therefore, we believe the best-fitting model is acceptable.

6. Discussion

Our estimates of the effective viscosities obtained using the late post-seismic deformation are compared with the effective viscosities reported in previous studies obtained using early post-seismic deformation covering different periods. Considering that the observed InSAR data are more sensitive to the lower crustal viscosity, the various lower crustal viscosity estimates are mainly compared. For the studies assuming a Burgers body, only the steady-state viscosity was used. For the studies considering the material difference between the northern and southern regions of the Kokoxili fault, the average value of the two sides was used.

For the Kokoxili earthquake (

Table 5. Estimated viscosities based on various observations following the Kokoxili earthquake.

	Earthquake	Viscoelastic Body	Post-Seismic Deformation Data	Viscosity (Pas)
1. Zhang [40]	Kokoxili	SLS	GPS data from 0–1 year after the earthquake	$\mathsf{\eta }=4\times {10}^{18}$
2. Diao [81]		Maxwell	GPS data from 0–4 months after the earthquake	${\mathsf{\eta }}_{\mathrm{lc}}\ge 1\times {10}^{18};{\mathsf{\eta }}_{\mathrm{um}}=1\times {10}^{18}$
3. Ryder [2]		Burgers	InSAR data from 2–5 years after the earthquake	$\mathsf{\eta }=4\times {10}^{19}$
4. Wen [41]		Maxwell	InSAR data from 2–6 years after the earthquake	${\mathsf{\eta }}_{\mathrm{lc}}=\left(2-5\right)\times {10}^{19};{\mathsf{\eta }}_{\mathrm{um}}=2\times {10}^{19}$
5. He [3]		Burgers	GPS data from 0–6.5 years after the earthquake	$\mathsf{\eta }=4\times {10}^{19}$
6. Shao * [42]		Maxwell	GPS data from 0–7 months after the earthquake	${\mathsf{\eta }}_{\mathrm{lc}}=9\times {10}^{18};{\mathsf{\eta }}_{\mathrm{um}}=1\times {10}^{20}$
7. This study		Maxwell	InSAR data from 14–19 years after the earthquake	${\mathsf{\eta }}_{\mathrm{lc}}=1.78\times {10}^{19};{\mathsf{\eta }}_{\mathrm{um}}=1.0\times {10}^{20}$

$\mathsf{\eta }$: viscosity for the viscoelastic half-space. *: the study considering the material difference between the northern and southern regions of the Kokoxili fault.

and Figure 7), the effective viscosity gradually increases with time during the post-seismic deformation process, consistent with the stress-dependent behavior inferred from laboratory experiments [81,82]. The effective viscosities obtained in this study are also similar to the steady-state viscosity calculated by Ryder [2] and Penchao [3], who used a Burgers body, which means that the viscoelastic layer under the Kokoxili earthquake has also almost reached its steady state.

7. Conclusions

This study calculated the late post-seismic deformation velocity of the 2001 Mw 7.8 Kokoxili earthquakes for the period between November 2014 and July 2021, which is much later than that reported in previous studies (generally covering the period until 6.5 years after these events), from Sentinel-1 A/B ascending and descending data. The deformation velocity field caused by the interseismic slip along the major faults in Tibet was first simulated. By comparing the simulated deformation velocity with the observed one, the maximum ratio of the simulated velocity to the observed velocity was found to be 42%, indicating continuity in the viscoelastic relaxation caused by the 2001 Kokoxili earthquake, and the interseismic signal cannot be ignored during the inversion.

By using the mixed model comprising the viscoelastic relaxation mechanism and the buried elastic dislocation model, the rheology of the northern Tibet Plateau lithosphere during this late post-seismic period was investigated and compared with previous studies that used early post-seismic deformation. The following conclusions can be drawn from this study:

• The best estimated viscosities for the lower crust and upper mantle were ${\mathsf{\eta }}_{\mathrm{lc}}=1_{-0.44}^{+0.78}\times {10}^{19} \mathrm{Pas}$ and ${\mathsf{\eta }}_{\mathrm{um}}=1_{+0}^{+0.78}\times {10}^{20} \mathrm{Pas}$, respectively. The estimated effective lower crustal viscosity was similar to the steady-state viscosity estimated by Ryder and Penchao, who considered a Burgers body. In addition, the temporal trend in the effective viscosity obtained by combining our results with those of previous studies was a significant logarithmic curve. Therefore, it can be concluded that the viscosity of the viscoelastic layer beneath the Kokoxili region has almost reached its steady state.
• The lower crustal steady-state viscosity of Songpan–Ganzi and Qiangtang terrane was ${10}^{19} \mathrm{Pas}$, which was larger than that obtained by the channelized lower crustal flow model (${10}^{17} \mathrm{Pas}$). This comparison means that the channelized lower crustal flow cannot be the main mechanism for the north Tibet Plateau deformation.

In addition, there are two aspects that require further study. The first is that the lateral heterogeneity in rheological structure between Songpan–Ganzi terrane and Qaidam Basin was not considered. The other is that a simplified, uniform, buried dislocation that represents the Kunlun fault was used because a precise interseismic aseismic slip rate was beyond the scope of the paper. Therefore, future work can be directed toward considering the lateral heterogeneity and using the multiple buried dislocations to invert the rheology structure of northern Tibet Plateau and precise interseismic slip rate along major faults in Tibet.