Rupture Kinematics and Coseismic Slip Model of the 2021 Mw 7.3 Maduo (China) Earthquake: Implications for the Seismic Hazard of the Kunlun Fault

Abstract

The 21 May 2021 Maduo earthquake was the largest event to occur on a secondary fault in the interior of the active Bayanhar block on the north-central Tibetan plateau in the last twenty years. A detailed kinematic study of the Maduo earthquake helps us to better understand the seismogenic environments of the secondary faults within the block, and its relationship with the block-bounding faults. In this study, firstly, SAR images are used to obtain the coseismic deformation fields. Secondly, we use a strain model-based method and steepest descent method (SDM) to resolve the three-dimensional displacement components and to invert the coseismic slip distribution constrained by coseismic displacement fields, respectively. The three-dimensional displacement fields reveal a dominant left-lateral strike-slip motion, local horizontal displacement variations and widely distributed near-fault subsidence/uplift deformation. We prefer a five-segment fault slip model, with well constrained fault geometry featuring different dip angles and striking, constrained by InSAR observations. The peak coseismic slip is estimated to be ~5 m near longitude 98.9°E at a depth of ~4–7 km. Overall, the distribution of the coseismic slip on the fault is highly correlated to the measured surface displacement offsets along the entire rupture. We observe the moderate shallow slip deficit and limited afterslip deformation following the Maduo earthquake, it may indicate the effects of off-fault deformation during the earthquake and stable interseismic creep on the fault. The occurrence of the Maduo earthquake on a subsidiary fault updates the importance and the traditional estimate of the seismic hazards for the Kunlun fault.

1. Introduction

Due to ongoing continental collision of the Eurasian and Indian plates, the Tibetan plateau is one of the regions globally with intense seismicity resulting from significant tectonic loading and crustal deformation. A series of large-scale active strike-slip faults and active sub-blocks are widely distributed across the Tibetan plateau. Among them, the Bayanhar block (Figure 1a), on the northern Tibetan plateau, is demonstrated to be the most seismically active sub-block in recent years. A series of major earthquake sequences occurred at the boundary of the Bayanhar block, including the 1997 Ms7.5 Manyi earthquake, the 2001 Ms8.1 Kokoxili earthquake, the 2008 Ms8.0 Wenchuan earthquake, the 2010 Ms7.1 Yushu earthquake, the 2013 Ms7.0 Lushan earthquake and the 2017 Ms7.0 Jiuzhaigou earthquake (Figure 1a). The seismic hazards of the existing unruptured sections of the main boundary faults, such as eastern part of the Kunlun fault, have been receiving extensive interest [1,2,3]. However, the seismic hazard of several important secondary faults in the interior of the block lacks detailed kinematic investigations and receives little attention due to limited geodetic observations. On 21 May 2021, a moment magnitude (Mw) 7.3 earthquake occurred in Maduo County (Qinghai Province, China). The epicentre was located at (34.59°N, 98.34°E), with a focal depth of 17 km, as determined by the China Seismological Network (http://www.csi.ac.cn/, accessed on 21 May 2021). This event demonstrates the great potential of generating major earthquakes of the subsidiary faults inside the blocks.

The Mw 7.3 Maduo earthquake occurred on a secondary fault within the eastern part of the Bayanhar block on the north-central Tibetan plateau [4]. Geometrically, the seismogenic fault of the Maduo earthquake is approximately subparallel to the main trace of the Kunlun fault to the north, near the big bend (97–99°E) of the middle-eastern section of the Kunlun fault (Figure 1b). According to the seismic records reported by the Chinese Earthquake Network Center, only a few historical earthquakes with M > 3 have occurred near Maduo in the last decade (Figure 1b), with the largest magnitude being Mw 5.2 in 2015 [4]. Additionally, the unanticipated Maduo earthquake was another major earthquake with M > 7 on the Tibetan plateau following the 2017 Jiuzhaigou earthquake. The seismic potential of the subsidiary fault in the Bayarhar block has been estimated to be low and not significant compared with the Kunlun fault to the north and the Ganzi-Yushu-Xianshuihe fault to the south [1,5,6]. It is still unclear whether the occurrence of the Maduo earthquake and the kinematics of the seismogenic fault are modulated by the activity of the Kunlun fault, or the stress accumulation and release of the subsidiary fault may influence the stress state of the main Kunlun fault.

Although the boundary of the Bayanhar block is a seismic zone capable of generating major earthquakes resulting from large fault slip rates (~6–11 mm/year) [6,7,8,9,10,11,12,13] and locking depths (>15 km), such as the Kunlun fault on the northern boundary, the seismic hazard of the large-scale subsidiary faults within the block has been overlooked and ignored to a great extent, which is further challenged by the difficulty of conducting field investigations on the higher plateau. The potential of major earthquakes within the block largely depends on the size of locked asperity on the fault, the interseismic fault slip rate, geometry and the stress perturbation induced by surrounding ruptures [14,15,16], which are important inputs of the seismic hazard model. The 2021 Mw 7.3 Maduo earthquake has been reported to have caused significant coseismic deformation and clear surface ruptures [17]. Geodetic measurements of coseismic deformation and the detailed analysis of fault ruptures caused by the Maduo earthquake may provide new insights into the fault kinematics of other active faults inside the block. Such information is also key to reveal the segmentation of the active faults within the Baryanhar block, the spatio-temporal heterogeneity of seismicity, and their relationship with the boundary faults.

In this study, we focus on determining the surface deformation characteristics and the geometry of the seismogenic fault of the 2021 Mw 7.3 Maduo earthquake. We use SAR images to map coseismic deformation fields and surface ruptures, to investigate the rupture kinematics and to constrain fault geometry parameters as well as coseismic slip distributions of the Maduo earthquake.

2. Data and Methods

2.1. InSAR Data and Processing

In order to measure the coseismic deformation produced by the Maduo earthquake, we apply conventional D-InSAR and offset-tracking methods to process Sentinel-1 SAR images, acquired in Terrain Observation with Progressive Scans (TOPS) mode from the European Space Agency (ESA), to obtain ascending and descending displacement fields in the line-of-sight (LOS) and ground range direction. The reference images were acquired on May 20, 2021 and the secondary images were acquired on May 26, 2021 for both ascending and descending tracks. See

Table 1. Parameters of the Sentinel-1A SAR images used in this study.

Orbit	Track	Acquisition Date		Perpendicular Baseline (m)	Incidence Angle (°)	Heading (°)
		Reference Image	Secondary Image
Ascending	106	20 May 2021	26 May 2021	−0.36	33~43	−13
Descending	99	20 May 2021	26 May 2021	3.22	33~43	−144

for other detailed parameters. All of the SAR images are interferometrically processed using the GAMMA software package [21], and multilook factors of two and ten are used in the azimuth and range directions, respectively, in order to filter phase noises and improve the signal-to-noise ratio (SNR) of the interferograms. Due to the significant along-track doppler centroid variation of TOPS SAR images, we use the enhanced spectral diversity (ESD) algorithm in overlapping areas across adjacent bursts to remove the residual azimuth phase ramp [22]. During interferogram formation, the Shuttle Radar Topography Mission (SRTM) 3-arc-second (90 m) DEM is used to remove the topographic phase. We filter and unwrap the interferograms using an adaptive filter [23] and minimum cost flow (MCF) [21] algorithm, respectively. The unwrapped interferograms are then geocoded to geographic coordinates.

We applied the intensity offset tracking algorithm, implemented in the GAMMA software, to derive the ground range displacement field using both ascending and descending images. We used a matching window of 200 by 40 pixels in range and azimuth directions, with an oversampling factor of 2 for improving the quality of the tracking results. The sampling steps are of 20 and 4 pixels in range and azimuth directions, respectively.

2.2. Three-Dimensional Displacement Field Decomposing

Due to the squint observation mode, the measurements of surface displacements using D-InSAR method are typically limited to the only one-dimensional LOS direction [24]. Theoretically, three or more InSAR displacement measurements with distinct perspectives, i.e., azimuth of the satellite and incidence angle of SAR images, are required to robustly resolve the three-dimensional components of surface deformation. However, this requirement is not always fulfilled in that data sources of SAR images are usually limited. To address such limitations, we use a combination of LOS and ground range observations of deformation on both ascending and descending tracks constrained by a strain model based algorithm [25] to reconstruct the three dimensional displacement fields of the Maduo earthquake. Particularly, this method is demonstrated to have the ability to considerably improve the accuracy of the solution in the north-south (NS) direction [24].

The basic strategy and formula are described as follows. We assume that a point of interest $P^0$, its coordinates and 3D components can be expressed as $p^0={\left[x^0{y}^0{z}^0\right]}^T$, $d^0={\left[e^0{n}^0{u}^0\right]}^T$, respectively. The coordinates and 3D components of the N points $P^i\left(i=1,2,\dots ,N\right)$ around it can be expressed as $p^i={\left[x^i{y}^i{z}^i\right]}^T$, $d^i={\left[e^i{n}^i{u}^i\right]}^T$, respectively. By using a strain model, the relationship between $d^i$ and $d^0$ can be expressed as:

$$d^i=H·\Delta p^i+d^0,i=1,2,\dots ,N$$

(1)

where $\Delta p^i=p^i-p^0={\left[\Delta x^i\Delta y^i\Delta z^i\right]}^T$, indicating the difference of coordinates between $p^i$ and $p^0$. H denotes the strain parameter matrix and can be divided into a strain tensor component ${\epsilon }_{ij}$ and a rotation tensor ${\omega }_{ij}$:

$$H={\epsilon }_{ij}+{\omega }_{ij}=\left[\begin{array}{ccc}{\epsilon }_{11}& {\epsilon }_{12}& {\epsilon }_{13}\\ {\epsilon }_{21}& {\epsilon }_{22}& {\epsilon }_{23}\\ {\epsilon }_{31}& {\epsilon }_{32}& {\epsilon }_{33}\end{array}\right]+\left[\begin{array}{ccc}0& {\omega }_{12}& {\omega }_{13}\\ -{\omega }_{12}& 0& {\omega }_{23}\\ -{\omega }_{13}& -{\omega }_{23}& 0\end{array}\right]$$

(2)

The strain tensor is symmetric, so only six components are independent. Symmetry requires that ${\epsilon }_{12}={\epsilon }_{21},{\epsilon }_{13}={\epsilon }_{31},{\epsilon }_{23}={\epsilon }_{32}$. In a simplified form, Equation (1) can be re-written as:

$$d^i=A_{sm}^{i}x,i=1,2,\dots ,N$$

(3)

where $A_{sm}^i$ is the strain model matrix, which is described in Equation (4). The matrix of unknown parameters $x={\left[e^0,n^0,u^0,{\epsilon }_{11},{\epsilon }_{12},{\epsilon }_{13},{\epsilon }_{22},{\epsilon }_{23},{\epsilon }_{33},{\omega }_{12},{\omega }_{13},{\omega }_{23}\right]}^{\mathrm{T}}$.

$$A_{sm}^i=\left[\begin{array}{cccccccccccc}1& 0& 0& \Delta x^i& \Delta y^i& \Delta z^i& 0& 0& 0& 0& \Delta z^i& -\Delta y^i\\ 0& 1& 0& 0& \Delta x^i& 0& \Delta y^i& \Delta z^i& 0& -\Delta z^i& 0& \Delta x^i\\ 0& 0& 1& 0& 0& \Delta x^i& 0& \Delta y^i& \Delta z^i& \Delta y^i& -\Delta x^i& 0\end{array}\right]$$

(4)

The relationship between the displacement field in LOS and the three dimensional displacement components at point $P^i$ can be mathematically expressed as:

$$L^i=A_{geo}^i·d^i$$

(5)

where $L^i=\left[\begin{array}{cccc}{L}_{LOS}^a& L_{LOS}^d& L_R^a& L_R^d\end{array}\right]$ is the LOS and ground range direction displacement of ascending and descending track of the point $P^i$, $A_{geo}^i$ is the transformation matrix between displacement fields and three dimensional components, which is shown in Equation (6):

$$A_{\mathrm{geo}}^i=\left[\begin{array}{l}\begin{array}{ccc}{a}_1^i& b_1^i& c_1^i\\ a_2^i& b_2^i& c_2^i\end{array}\\ \begin{array}{ccc}{a}_3^i& b_3^i& c_3^i\\ a_4^i& b_4^i& c_4^i\end{array}\end{array}\right]$$

(6)

$$\{\begin{array}{l}{a}_1^i=-\sin \left({\theta }_{as}^i\right)\sin \left({\alpha }_{as}^i-3\pi /2\right),b_1^i=-\sin \left({\theta }_{as}^i\right)\cos \left({\alpha }_{as}^i-3\pi /2\right)\\ a_2^i=-\sin \left({\theta }_{des}^i\right)\sin \left({\alpha }_{des}^i-3\pi /2\right),b_2^i=-\sin \left({\theta }_{des}^i\right)\cos \left({\alpha }_{des}^i-3\pi /2\right)\\ a_3^i=\sin \left({\alpha }_{as}^i-3\pi /2\right),b_3^i=\cos \left({\alpha }_{as}^i-3\pi /2\right)\\ a_4^i=\sin \left({\alpha }_{des}^i-3\pi /2\right),b_4^i=\cos \left({\alpha }_{des}^i-3\pi /2\right)\\ c_1^i=\cos \left({\theta }_{as}^i\right),c_2^i=\cos \left({\theta }_{des}^i\right),c_3^i=c_4^i=0\end{array}$$

(7)

where $\alpha$ is the heading angle of satellite, and $\theta$ is the incidence angle for a specific pixel.

Equations (3) and (5) can be combined as:

$$\begin{array}{c}{L}^i=A_{geo}^i·d^i=A_{geo}^i·A_{sm}^{i}x,i=1,2,\dots ,N\\ \underset{4K\times 1}{L}=\underset{4K\times 12}{A}·\underset{12\times 1}{x}\end{array}$$

(8)

2.3. Slip Distribution Inversion Method

We use the widely acknowledged steepest descent method (SDM) [26] to invert for the coseismic slip distribution of the Maduo earthquake to investigate the relationship between kinematics of surface ruptures and deep slip on the fault. The basic mathematical description of our inversion can be summarized as follows:

$$f(m)={\displaystyle \sum _{i=0}^n\Vert Gm(}x)-d{\Vert }^2+{\alpha }^2\Vert Hm{\Vert }^2\to \min$$

(9)

where G is the Green’s function for the elastic half-space, m is the matrix of the slip vector, d is the matrix of geodetic measurements on the surface, ɑ is the smoothing factor, H is the Laplacian operator on slip, and $\Vert Hm{\Vert }^2$ is the slip roughness.

A grid search method [27] is adopted in our inversion to find the best fitting fault geometry of the seismogenic fault, due to the fact that the fault geometry in our model could be simplified and described by a few parameters, including the bottom/top depths of fault segments and the dip angle. Owing to the fact that the top and bottom depths of the fault are fixed at 0 and 20 km, respectively, in our inversion, the only free parameter related to the fault geometry is the dip angle of fault segments. We assume a uniform dip angle for all the depths in that no evidence and observations can illustrate a curved fault geometry.

In our inversion, the fault plane is discretized into a series of 2 km × 2 km sub-fault patches. We use both ascending and descending InSAR measurements to constrain the fault geometry, because of their sensitivity to local fault geometries. To decrease the number of points from the InSAR observations and to maximize the computational efficiency, we employ the uniform sampling method to downsample ascending and descending displacement fields to a manageable amount. We calculate the Green’s functions for all the dislocations embedded in a homogeneous elastic half-space using a Poisson’s ratio of 0.25. The fault surface traces are geodetically extracted based on obvious displacement discontinuities, which is nicely consistent in both ascending and descending interferograms, and we divided the fault into 5 segments. We assume a uniform dip angle of five fault planes from geodetically detected fault traces. Informed by the preliminary analysis of seismic observations, we allow the rake angle of the coseismic slip to vary in the range of −45°~45° for all the fault segments in the model, which is corresponding to the strike-slip motion with possible thrust or normal slip on the fault. The fault plane extends to a maximum depth of 20 km, because preliminary inversion indicates that the depth of 20 km is sufficient to resolve all the coseismic slip during this event, and most of the coseismic slip is distributed in the depth of shallower than 20 km. After fixing the top and bottom depth, we only need to search for five dipping angles, which makes the grid search method practicable. We test a suite of fault models with a series of dip angles, invert the corresponding slip distributions, and compute associated root-mean-square between data and model predictions. The Root Mean Square (RMS) misfit is calculated using the following expression:

$$RMS=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m{\left(d_{obs}^{ij}-d_{\mathrm{mod}}^{ij}\right)}^2}}}{nm}}$$

(10)

where $d_{obs}^{ij}$ and $d_{\mathrm{mod}}^{ij}$ are the observed and modeled surface displacements, respectively.

3. Results

3.1. Coseismic Deformation Fields

Figure 2a,c show our derived coseismic interferograms on the ascending and descending tracks. Overall, InSAR observations captured the complete and large-scale coseismic deformation field of the Maduo earthquake. The nearly linear and continuous displacement discontinuity only exists in the near-fault area resulting from the low coherence produced by the large displacement gradient across surface ruptures, which informs us the possible trace of surface ruptures. Additionally, InSAR observations, on both ascending and descending tracks, demonstrate a slightly asymmetrical pattern across the fault, indicating a possible fault geometrical variation along the whole rupture. In our displacement field, a positive displacement value indicates surface motion towards the satellite, while a negative value means surface motion away from the satellite. The northern part of the deformation field shows negative displacement in the descending deformation maps (Figure 2b), while positive displacements in the ascending deformation maps (Figure 2d). Such result demonstrates a left-lateral strike-slip motion of the seismogenic fault, which is generally in agreement with focal mechanism solutions (

Table 2. Focal mechanism parameters of the Maduo earthquake.

Source	Mw	Focal Depth (km)	Longitude (°)	Latitude (°)	Strike (°)	Dip (°)	Rake (°)
GCMT	7.4	12	98.46	34.65	282	83	−9
USGS	7.3	10	98.2458	34.6125	92	67	−40

). InSAR deformation maps clearly show that the whole displacement field covers an area of approximately 200 km × 70 km. The displacement fields on both the ascending and descending tracks indicate a maximum peak-to-trough offset of ~2 m in the LOS direction near the central segment of the fault (Figure 2g,i).

Figure 2e,f show the coseismic displacement fields in the ground range direction resolved by the offset tracking method, and the incoherent regions to the west of the fault are water bodies. The maximum peak-to-trough displacement on the ascending track is estimated to be ~1.2 m (Figure 2h), while the maximum peak-to-trough displacement on the descending track is about ~1.3 m (Figure 2j).

All the InSAR observations with high quality reveal a generally consistent and rather clear surface rupture trace with an WNW orientation and a length of more than ~160 km. We find obvious geometric variations near both the east and west end of the surface rupture, which may control the termination of dynamic rupture during the Maduo earthquake. Figure 2g–j show coseismic displacement distributions along the rupture zone (AA′ profile) in the LOS and ground range direction, respectively. To first order, it suggests that the coseismic displacement in the near-field of the fault is heterogeneous and multi-peaked (Figure 2g,i), which is the product of the complexity of fault geometry (i.e., step-overs, fault bends), fault segmentation along strike and inelastic off-fault deformation. InSAR captured displacement pattern along strike is also consistent with field investigations immediately conducted after the earthquake.

3.2. Complete Three Dimensional Displacement Fields

Our derived three-dimensional displacement fields of the Maduo earthquake are illustrated in Figure 3. Our result confirms that the Maduo earthquake is dominated by left-lateral strike-slip motion, suggested by the dominant contributions from horizontal, i.e., east-west (EW) and north-south (NS), displacement components. The amplitude (exceeding ±0.8 m) of the EW displacement is the largest, and the EW displacement field is largely associated with the left-lateral strike-slip shear of the seismogenic fault. The EW displacement field also contributes to the observed large displacement gradient near the fault (Figure 3a). The NS component has a maximum displacement of ~0.8 m in the northern side of the fault and ~0.9 m in the southern side of the fault. The spatial distribution of the NS displacement field is more complex and spatially heterogeneous with some small-scale near-fault features at different fault segments, which are spatially corresponding to the local fault geometry (Figure 3b). The vertical displacement component has the smallest amplitude (<0.5 m), only about half of the horizontal components (Figure 3c). Moreover, the vertical displacement distribution is very uneven.

The horizontal displacement vector field depicts the overall displacement seismogenic fault movement pattern, direction and magnitude of both the near- and far-field deformation of the Maduo earthquake. It indicates that the northern and southern sides of the fault as a whole is rotating clockwise (Figure 3c), which is consistent with the general pattern of coseismic deformation of a long strike-slip fault. We also find a slightly asymmetrical displacement pattern with a larger displacement magnitude in the northern side of the fault, which is also consistent with the widely distributed uplift in the northern side of the seismogenic fault. Such features are uniquely related to the geometry of the north-dipping fault, which is confirmed by our following geodetically constrained fault slip model in Section 3.3.

3.3. Fault Geometry and Coseismic Slip Distributions

RMS versus dip angle trade-off curve can be used (Figure 4a–e) to locate the best-fitting dip angle which best explains the observations. The optimal smoothing factor is selected based on the trade-off between the misfit and model roughness (Figure 4f). Our result indicates that the seismogenic fault of the Maudo earthquake is highly segmented, as revealed by relatively large variations of the dip angle along the whole length of the fault. Overall, the fault is north dipping with dip angles of 70°, 65°, 80°, 80°, 85° for fault segments 1–5, respectively. The north-dipping fault geometry is also consistent with the clear offset between the rupture trace on the surface and the distribution of the deep relocated aftershocks (Figure 1b and Figure 3c) and the widely distributed uplift deformation to the north of the fault (Figure 3c). East of the triple junction (segments 4 and 5), our fault model somewhat conflicts with observations: segment 4 dipping to the north, segment 5 dipping to the south, and vertical displacement showing slight subsidence in this area. Since faults are sub-vertical here, it may be challenging to derive the dip direction of fault segment 5. For other fault segments, we do not particularly emphasize the respective dipping direction of the fault, but still referred to them as sub-vertical faults, which is consistent with the dominated strike-slip motion during the Maduo earthquake, and the slightly asymmetric displacement fields observed by InSAR data (Figure 2).

Our preferred coseismic slip distributions of the Maduo earthquake are demonstrated in Figure 5. Our inversion result constrained by ascending and descending InSAR measurements suggests relatively shallow rupture and several large asperities along five different fault segments (Figure 5). Overall, the coseismic slip during this event is generally confined in the depth range of 0–15 km, indicating that the coseismic rupture did not substantially penetrated into the shear zone or the much deeper part of the seismogenic depth (>15–20 km). We find the most significant left-lateral slip is distributed near the triple junction area of the fault (~99.2°E) with a peak magnitude of ~5 m at a depth of 6–10 km. The total seismic moment is ~1.45 × 10²⁰ N m, which is equivalent to an Mw 7.41 earthquake, consistent with the focal mechanism released by Global CMT catalog. Our preferred fault model gives a satisfied fit to the ascending and descending InSAR data (Figure 6), which can nicely explain the downsampled InSAR data with RMS values of 8.8 mm and 6.7 mm for the ascending and descending InSAR data, respectively. Large residuals only exist in the near-field of the fault, which possibly originates from the inelastic off-fault deformation, over-simplification of fault geometry near the surface and some other phase noises associated with incoherence in the interferograms.

To investigate the relationship between fault geometry and rupture kinematics, we calculate the coseismic displacement offsets along the surface rupture using our derived three-dimensional displacement fields (Figure 3). The displacement offsets are calculated by the difference between the mean value of several points (ranging from 30 to 100 pixels with a resolution of ~30 m) on both sides of the surface ruptures. Although the observations here do not explicitly reflect the real coseismic offsets across the ruptures in nature, limited by the resolutions of InSAR observations and smoothness during the calculation, the investigations mainly aim to understand the first-order feature of both the deep, shallow slip on the fault and the deformation on the surface. We quantify fault geometry using the geometrical curvature and fault strike. These two types of parameters are useful to indicate the small- and large-scale fault geometry variations, which may control short- and long-wavelength surface deformation. We measure the coseismic offsets along the fault rupture in north-south, east-west and fault-strike parallel directions. We also compare surface rupture kinematics with our inverted coseismic slip on the fault plane. All the results are illustrated in Figure 7.

To first order, the distribution of the coseismic offsets in the north-south direction (Figure 7c) is much more heterogeneous than that in the east-west displacement (Figure 7b). The rough distribution of the north-south displacements is remarkably correlated to the local fault strike change as inferred from the fault curvature (Figure 7a). Local fault strike change, leading to the fault curvature variation, results from step-overs, fault bends and flower structures of the strike-slip fault. These tectonic structures dominantly control the near-fault and on-fault deformation and may change the local deformation pattern substantially. Our InSAR measured coseismic offsets cannot directly be compared with offset measurements from field investigations, due to the fact that the offset measurements only capture on-fault displacements while InSAR data measure both the on-fault and the mid- to far-field coseismic offsets [28]. The coseismic displacement offset along fault strike (Figure 7d) has similar patterns but larger magnitude than coseismic offset in the east-west direction (Figure 7b), and it is expected that the east-west component dominantly contributes to the measured coseismic offset along the fault strike. The comparison between coseismic offset on the surface (Figure 7d) with inverted coseismic slip on the fault plane (Figure 7e) suggests that deep slip and shallow fault offsets are also highly correlated, and greater deep slip generates larger surface offsets (>0.5 m). We also find that the distribution of aftershocks is spatially complementary with the coseismic slip, with aftershock clusters surrounding multiple asperities.

4. Discussion

4.1. Depth Distribution of Coseismic Slip and Aftershocks

One of the interesting kinematic features of the Maudo earthquake is the long (~160 km) and relatively shallow rupture and the fact most resolvable coseismic slip is confined at the depth of shallower than 13–15 km (Figure 5 and Figure 7). We summarize the depth distribution (moment profiles) of coseismic slip and aftershocks along five segments in Figure 8. It indicates a complementary distribution between relatively shallow coseismic slip and deep triggered aftershocks. Although the coseismic rupture breaks to the surface, the obvious but not significant shallow slip deficit is still observed for all five fault segments. The shallow slip deficit during the earthquake likely stems from distinct friction properties (velocity weakening and velocity strengthening) at various depths of the seismogenic zone. Under this hypothesis, the accumulated strain near the surface during the interseismic period may be released via aseismic stable creep and/or afterslip following the earthquake, or more distributed off-fault deformation triggered by coseismic rupture [29,30,31,32,33,34].

It is generally accepted that afterslip and aftershocks are favored or triggered by the static stress change imposed by the main shock [34,35,36]. To assess the effect of coseismic slip and associated stress changes on the occurrence and spatial distribution of afterslip and aftershocks, we calculate the change of Coulomb Failure Stress (CFS) based on our preferred coseismic slip model assuming an effective friction coefficient of 0.4. The result is shown in Figure 5c. Overall, the correlation between the positive CFS change and aftershock distribution is strong on all five fault planes. This suggests that the Coulomb stress change may play an important role in controlling the spatial distribution and evolution of aftershocks in this case.

Here, we also obtain the postseismic interferograms using SAR images acquired on 26 May 2021 (six days after the earthquake) and 1 June 2021 (twelve days after the earthquake), which captures the early short-term postseismic deformation (Figure 9). We find that the detectable postseismic deformation only occurred along the section between longitude ~98°and ~98.6° (~50 km), and only moderate coseismic slip (~2–3 m between 0 and 10 km) occurred along this segment (Figure 7e). The scale (near-field), magnitude (<5 mm accumulated in six days), wavelength and spatial pattern (also left-lateral strike-slip motion) suggest that the postseismic deformation is likely produced by shallow afterslip, which have the potential to partially compensate the shallow slip deficit during the coseismic rupture. Considering that this paper mainly focuses on coseismic studies, we do not intend to simulate post-earthquake deformation in this study. However, the limited postseismic deformation compared with long coseismic rupture (~160 km) and significant coseismic stress perturbation below the depth of 13–15 km (>2 Mpa, Figure 5c) arises further questions related to the mechanism and process of stress relaxation after the Maduo earthquake. In the follow-up study, we will use more intensive long time post-seismic InSAR data to carry out a detailed study.

4.2. Implications on Seismic Hazard Estimate for the Kunlun Fault

Traditionally, considerable efforts on the estimate and assessment of the seismic hazard along the Kunlun fault focus on the Xidatan-Dongdatan (~93°–95°) and Maqin-Maqu (~99°–103°) segments, owing to the well determined fast fault slip rate (~6–12 mm/year) and nearly fully locked status, indicating a possible phase at the late earthquake cycle. Additionally, the almost constant distribution of fault slip rate along the Kunlun fault between 92°and 100° provides evidence for the inference of the block-like behavior of the Baryanhar block, even under an inhomogeneous tectonic loading along the entire Kunlun fault, potentially compensated for by a heterogeneous frictional resistance [10,11]. These results are compatible with the view that the seismic hazard of unruptured fault segments on the Kunlun fault is high and significant, while the earthquake potential inside the Bayanhar block is low. However, the occurrence of the Maduo earthquake challenges the widely accepted view and the Bayanhar block may deform in a more distributed way through the strain accumulation and release on a series of sub-parallel and large-scale strike-slip faults (Figure 1). Due to the spatial proximity between these faults, the block models used to image the elastic coupling and large locked asperities on the fault need to address the potential trade-off of the fault slip rates and elastic coupling distributions between the Kunlun fault and these subsidiary faults.

Another important implication is related to the earthquake cycle modelling for the Kunlun fault. Contemporaneous slip-rates (decadal scale) are typically established with geodetic observations using a series of pre-defined faults and block models [12,14]. Long-term fault slip rates on the Kunlun fault is also well determined by geological survey [10,11]. However, near the big bend of the Kunlun fault, these two types of estimates have discrepancy, which has been tentatively explained by earthquake cycle effects [19]. It is important to understand the inconsistency between the geodetic and geological slip rates on various segments of the Kunlun fault. The extensive debates on the inconsistency between two types of measurements are associated with the Tuosuo Lake segment of the Kunlun fault [19], which is the largest releasing bend of the Kunlun fault, also near the seismic zone of the Maduo earthquake. The incorporation of more complex earthquake cycle models, such as the model incorporating a clustered earthquake sequence, are also highly required.

5. Conclusions

The 2021 Mw 7.3 Maduo earthquake occurred in the interior of the Baryanhar block near the big bend of the Kunlun fault. The occurrence of the Maduo earthquake challenges the common view that the Bayanhar block largely behaves like a block as inferred from the generally constant fault slip rate distribution of the Kunlun fault along sections between longitude 92°E and 98°E. In this study, we map significant coseismic deformation using ascending and descending D-InSAR and offset-tracking observations. Multi-perspective observations facilitate reconstruction of the three-dimensional displacement fields. High-quality InSAR observations also help to map the clear surface trace of the coseismic ruptures. We use the geodetic measurements on the surface to investigate the fault geometry and rupture kinematics of the seismogenic fault to understand their relationship. Three-dimensional displacement field confirms that the Maduo earthquake is dominated by the left-lateral strike slip with the largest horizontal offset of up to ~1.5 m, which is generally consistent with field investigations [17]. The widely distributed uplift and subsidence deformation in the near-field of the fault is largely associated with fault geometry variations. Our preferred fault model indicates that the Maduo earthquake have the peak slip ~5.0 m at the depth of ~8 km near the triple junction of the eastern end of the fault. The event was characterized by dominantly left-lateral strike-slip motion with a slight normal dip-slip component. The relatively complex fault geometry is quite different from the boundary fault of the Baryanhar block, i.e., the long and sub-vertical Kunlun fault to the north. Some other large-scale subsidiary faults need more attention related to their seismic hazard. However, the fault slip rate and locking depth of these faults are difficult to constrain and quantified using interseismic geodetic observations because of more distributed deformation inside the block as well as long earthquake intervals.