# Physics-based assessment of earthquake potential on the Anninghe-Zemuhe fault system in southwestern China

Faqi Diao, Huihui Weng, Jean-Paul Ampuero, Zhigang Shao, Rongjiang Wang, Feng Long, Xiong Xiong

► To cite this version:

Faqi Diao, Huihui Weng, Jean-Paul Ampuero, Zhigang Shao, Rongjiang Wang, et al., Physics-based assessment of earthquake potential on the Anninghe-Zemuhe fault system in southwestern China. Nature Communications, 2024, 15 (1), pp.6908. 10.1038/s41467-024-51313-w. hal-04885055

HAL Id: hal-04885055

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

Submitted on 15 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Physics-based assessment of earthquake potential on the Anninghe-Zemuhe fault system in southwestern China

Received: 28 November 2023

Accepted: 5 August 2024

Published online: 12 August 2024

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

Check for updates

Faqi Diao Ⓔ<fnref n="1" />✉, Huihui Weng Ⓔ²✉, Jean-Paul Ampuero Ⓔ³, Zhigang Shao⁴, Rongjiang Wang<fnref n="1" />,⁵, Feng Long Ⓒ⁶ & Xiong Xiong<fnref n="1" />

The seismic hazard of a fault system is controlled by the maximum possible earthquake magnitude it can host. However, existing methods to estimate maximum magnitudes can result in large uncertainties or ignore their temporal evolution. Here, we show how the maximum possible earthquake magnitude of a fault system can be assessed by combining high-resolution fault coupling maps with a physics-based model from three-dimensional dynamic fracture mechanics confirmed by dynamic rupture simulations. We demonstrate the method on the Anninghe-Zemuhe fault system in southwestern China, where dense near-fault geodetic data has been acquired. Our results show that this fault system currently has the potential to generate Mw7.0 earthquakes with maximum magnitudes increasing to Mw7.3 by 2200. These results are supported by the observed rupture extents and recurrence times of historical earthquakes and the b values of current seismicity. Our work provides a practical way to assess the earthquake potential of natural faults.

Large earthquakes reoccur on natural faults due to the interaction between tectonic loadings and frictional resistance. The characteristic earthquake model on an isolated seismogenic zone is often used to conceptualize observations of recurring earthquakes whose magnitudes are similar, such as the moderate earthquakes in the Parkfield segment of the San Andreas fault<fnref n="1" />,² and small repeating earthquakes³. For large subduction earthquakes, the supercycle model is required to explain the more complicated rupture behaviors on faults where multiple large earthquakes of different sizes occur in one complete supercycle⁴. Predicting the recurrence times of large earthquakes is challenging, even for the simple characteristic model⁵,⁶, mainly due to the unpredictable nucleation conditions⁷ and the uncertainties in stress states⁸. While we still cannot predict where and when earthquakes will start, recent advances⁹,¹⁰ can help us pinpoint where earthquakes might stop, thus how large they might grow.

An accessible way to assess future seismic hazard is to estimate the maximum possible earthquake magnitudes on faults as a function of elapsed time, capturing for example the possibility that a fault might not be ready to break over its whole length. Kinematically-based models¹¹⁻¹³ or dynamically-based⁸,¹⁴⁻²⁰ models on various fault systems, including estimates of their slip budget, have been used to infer their potential seismic moments. Specifically, fault geometric complexities²¹,²², such as fault bends and step-overs, have been integrated into single rupture simulations¹⁵⁻¹⁸ and earthquake cycle models¹⁹,²⁰. Furthermore, complex fault systems have also been incorporated with a physics-based statistical simulator²³ to reproduce seismic hazard statistics in Southern California. Both kinematically-based and dynamically-based models provide important insights on rupture segmentation and seismic hazard assessment. However, kinematically-based methods, although efficient, may overestimate the maximum magnitudes of potential earthquakes due to the lack of

<footnote n="1">¹Hubei Subsurface Multi-Scale Imaging Key Laboratory, School of Geophysics and Geomatics, China University of Geosciences, Wuhan, China. ²State Key Laboratory for Mineral Deposits Research, School of Earth Sciences and Engineering, Nanjing University, Nanjing, China. ³Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Géoazur, Valbonne, France. ⁴Institute of Earthquake Forecasting, China Earthquake Administration, Beijing, China. ⁵GFZ German Research Centre for Geosciences, Potsdam, Germany. ⁶Seismological Bureau of Sichuan Province, Chengdu, China. ✉e-mail: fqdiao@cug.edu.cn; weng@nju.edu.cn</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Tectonic setting and fault coupling ratio on the Anninghe-Zemuhe fault (AZF) system. A Tectonic setting of the study area. Black solid lines indicate active faults obtained from the public Data catalog of the China Seismic Experiment Site. Red and gray triangles represent the 24 newly-developed near-field GNSS stations and existing stations used for constraining the kinematic coupling model, respectively. Cyan stars indicate historical events (M > 6.5) in the past 500 years (ref. 24). Blue vectors show the GNSS velocity in this region ("Data availability" section). ANHF: Anninghe fault; ZMHF: Zemuhe fault; DLSF: Daliang Shan fault; LJ-XJHF: Lijiang-Xiaojinhe fault; LMSF: Longmen Shan fault. Cyan vectors in the upper left inset show the block motions. Red star and ellipse represent the epicenter</figcaption></figure>

and rupture extent of the 2022 Mw 6.6 Luding earthquake. B Coupling ratio on the AZF and DLSF. The fault changes from a completely locked state to a fully steady slip state as the coupling ratio varies from 1 to 0. Gray dots show the seismicity distribution that was projected onto the faults. AZF-1 to AZF-6 show the identified asperities on the AZF. Purple dots represent the near-fault cities. The topography data used in the figure is obtained from NOAA Data Catalog. The GNSS velocities on the new stations are available at Table S1. The seismic catalog data is available in Supplementary Data 1. The figure was drawn using the GMT5 software ("Code availability" section).

integration with rupture dynamics, while dynamically-based numerical simulations are too computationally expensive for efficient hazard assessment and thorough sensitivity analysis.

The Xianshuihe-Xiaojiang fault (XXF) is one of the longest left-lateral faults in southwest China with a spatial scale (-1100 km) comparable to the San Andreas fault in California. The deformation in the middle part of the XXF is accommodated by two branches: the mature Anninghe-Zemuhe fault (AZF) and immature Daliangshan fault (DLSF) (Fig. 1). Paleoseismic evidence shows that large earthquakes frequently occurred on the AZF in the past 500 years, including six M > 6.5 events and three M > 7.5 events²⁴, whereas there is barely paleoseismic evidence of large events on the DLSF in the past 500 years²⁵. The northern segment of the AZF is now thought to be close to its reoccurrence period since the last M-7.5 event in 1480 (ref. 24). However, questions persist about the elastic energy currently accumulated on the two faults and about the maximum possible earthquake magnitudes the faults can accommodate currently and in the future.

Here, we combine kinematic inversions and a physics-based theoretical model based on three-dimensional (3D) dynamic fracture mechanics to show that the demonstrated fault system has already the potential to produce Mw > 7 earthquakes and the maximum possible magnitude will increase to Mw > 7.3 in 2200.

## Results

### Separated asperities on the Anninghe-Zemuhe fault system

Fault coupling controls the accumulation of slip deficit on faults and thus provides key constraints on the location and rupture extent of potential large earthquakes²⁶⁻²⁸. To obtain a high-resolution fault coupling distribution, we combine observations on 24 newly-developed near-field GNSS (Global Navigation Satellite System) stations and on previous stations around the AZF and DLSF (Fig. 1A) to estimate the long-term crustal velocities (Fig. 1 and Table S1) for the

kinematic inversion of fault coupling ("Methods" section). Because the AZF has a near N-S trending that limits the applicability of the InSAR (Interferometric Synthetic Aperture Radar) technique, the near-field GNSS observations in this study can provide valuable constraints on fault coupling distribution. In the interseismic stage, viscoelastic relaxation at depth can affect crustal deformations in middle- and far-filed areas²⁹. Ignoring this effect in the kinematic inversions results in larger slip deficits on deep fault segments to compensate, which thus leads to an overestimation of their coupling ratio³⁰. We therefore adopt viscoelastic deformation models to infer the coupling ratio of the AZF and the DLSF. Our results further confirm that the elastic model overestimates the fault coupling on deep patches (Fig. S1). We find that the inferred geodetic slip rate on the AZF decreases from north to south, from 6.1 mm/a to 4.9 mm/a, which agrees well with geologic observations of 6.5 ± 1.0 mm/a³¹. The DLSF has a slip rate of 3.3–4.8 mm/a, responsible for ~40% of the sinistral motion across the fault system, which is also consistent with the geologic estimation of ~3–4 mm/a³². The fault coupling ratio on the AZF shows "string-beads-shaped" lateral variation along the strike, featuring six high-coupling asperities with spatial scales of 20–40 km (Fig. 1B). Along the DLSF with a comparable length to the AZF, we identify four high-coupling asperities. Checkerboard tests show that asperities with scales >20 km can be effectively resolved (Fig. S2), owing to the newly developed dense near-field stations (Fig. 1A), which therefore provides sufficient resolution in our inverted coupling model.

The AZF and the DLSF have different fault maturities³²: the mature AZF has more continuous surface traces than the immature DLSF (Fig. 1A). Seismological³³,³⁴ and experimental³⁴ studies show that the fault topographic roughness, representative of the fault structural maturity, may control earthquake sizes: smoother fault interfaces tend to produce larger earthquakes. Here, our results show that the AZF and the DLSF have distinct patterns of coupling ratio distributions that may
also be caused by different fault maturities: the mature AZF has shorter
scales of coupling heterogeneities than the immature DLSF. In addi-
tion, lacking paleoseismic evidence of large earthquakes on the
immature DLSF²⁵ is qualitatively consistent with previous studies of
fault roughness³³,³⁵. If the experimental results of fault roughness in the
laboratory (e.g., ref. 35) can be upscaled to natural faults, the AZF may
tend to produce larger events than the DLSF, yet a conclusive com-
parison should also consider their accumulated slips, geological ages,
tectonic loading conditions, and other relevant factors.

Spatial correlation with historical ruptures and current seismic
b values

The laterally variable fault coupling is strongly related with the accu-
mulated shear stress and rupture extents of historical earthquakes. We
first can use the rupture extents of large historical earthquakes to
validate the reliability of the obtained fault coupling maps. As shown in
Fig. 2A, there is a strong spatial correlation between the inverted fault
asperities on the AZF and the observed rupture extents of large his-
torical earthquakes in the past -500 years, delineated by comprehen-
sive analysis of multi-disciplinary observations²⁴. Three M 6.8
earthquakes, which occurred in 1952, 1489, and 1732, only ruptured
single isolated asperities, whereas three M 7.5 events in 1480, 1536, and

1850 connected two adjacent asperities. If the rupture is confined within an isolated asperity, the earthquake magnitude is generally less than Mw 7, but it can reach up to Mw 7.5 if ruptures break through multiple adjacent asperities. Besides, we observe two low coupling zones, i.e., north of asperity AZF-3 and south of asperity AZF-4, that seem to be long-standing barriers as no ruptures breaking through them have been observed. The phenomenon that historical ruptures were generally impeded by low coupling areas strongly demonstrates that fault coupling controls the rupture extent and magnitude of large events. The correlation also indicates that the fault coupling may represent a long-term stable feature and thus we propose that our results based on decades of observations can be extended to at least one earthquake cycle.

We further calculate the stressing rate distribution based on the
fault coupling ratio on the AZF and estimate the current accumulated
shear stress on each asperity since the last large earthquakes (Fig. 2B;
“Methods” section). The *b* value has been suggested to act as a stress
meter that is inversely correlated with the stress level³⁶,³⁷. To compare
with the calculated stress, we derive the *b* value in this region based on
the updated earthquake catalog in the period from 2009 to 2023 and
draw it along the AZF (“Methods” section). As shown in Fig. 2C, we
observe two areas of low *b* values near asperity AZF-1 and AZF-4,

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

**Fig. 2 | Spatial variation of coupling ratio, historical rupture, b value, and stress along the Anninghe-Zemuhe fault (AZF).** A Comparison between the depth-averaged (0–15 km) coupling ratio (blue dots with gray error bars) and the rupture extent of large historical earthquakes (orange lines, ref. ²⁴) on the AZF. The gray error bars show 1-σ uncertainty of the coupling ratio inferred from statistical analysis (“Method” section). The solid orange lines indicate historical ruptures of high

reliability, whereas the dashed lines represent observations with relatively large
uncertainty. AZF-1 to AZF-6 at the top of the panel show the identified asperities on
the AZF. B Distribution of accumulated shear stress since 1480. C The *b* values along
the AZF, the error bar show 1-σ uncertainty. D Comparison between recurrence
intervals calculated based on the physics-based rupture theory (black line) and
those inferred from geologic observations (red bars).
<figure><img src="image_6.png" /><figcaption>**Fig. 3 | The maximum possible magnitude of earthquakes on the Anninghe-Zemuhe fault (AZF).** A The estimated distribution of the potential energy release rate and fracture energy along strike before the 1952 M 6.8 earthquake. The top histogram shows the seismicity distribution along the strike between 2009 and 2023. The colored band accounts for uncertainties of the $G_c$ scaling. B The maximum possible earthquake magnitude as a function of hypocenter location estimated by the 3D theory (purple symbols; vertical error bars indicate the</figcaption></figure>

uncertainties caused by the fitting of the $G_c$ scaling) and by numerical simulations (blue curve) with the same $G_c$ scaling. The thick red line marks the ruptured region of the 1952 M 6.8 earthquake. C, D Same as (A, B) but in 2023. In B and D the pink areas indicate the high-hazard regions where the maximum possible earthquake magnitudes are close to Mw 7 and the inset in (D) shows the maximum possible earthquake magnitude as a function of time until the year 2250.

indicating highly stressed conditions, which are consistent with the high shear stress on the fault accumulated since the last large events (Fig. 2B). On asperity AZF-3 and AZF-6, the *b* value is relatively high, revealing a relatively low-stress level due to short elapsed times on these asperities (Fig. 2B). These results highlight that the shear stress constrained by the geodetic and paleoseismic observations is qualitatively consistent with the *b* value based on the statistical properties of current seismicity.

## Maximum possible earthquake magnitudes

The accumulated slip deficit and shear stress, combined with 3D dynamic rupture theory, allow us to assess the maximum possible earthquake magnitude on the AZF and the DLSF. In dynamic rupture theory⁹, rupture propagation and arrest on a fault are controlled by the accumulated elastic energy $G_0$ from tectonic loadings and the dissipation energy $G_c$ due to frictional resistance. This theory predicts that rupture accelerates when $G_0 > G_c$ and decelerates when $G_0 < G_c$, encapsulating a physics-based arrest criterion for dynamic ruptures: $\int_{x_0}^{x_0+L} (1 - G_c(x)/G_0(x))dx = 0$, where *L* is predicted along-strike arrest location distance from the hypocenter (located at $x_0$). The criterion applies to each of the rupture fronts on bilateral ruptures. For the AZF (Fig. 3) and DLSF (Fig. S3), the distribution of $G_0(x)$ along strike can be derived based on the inverted coupling map and the elapsed time since the last large earthquake and $G_c(x)$ can be determined using a power-law slip-weakening friction model to fit global earthquake observations (“Methods” section).

We first use this theory to estimate the maximum possible earthquake magnitude on the AZF in 1952 before the occurrence of the M 6.8 earthquake, by assuming various hypocentre locations (Fig. 3B). We find that before the 1952 event, $G_0$ on the asperity AZF-3 is higher than $G_c$ (Fig. 3A), and thus it is in a favorable state for ruptures. The theoretical model shows that if a rupture nucleates around asperity AZF-3 the fault has a potential to generate a Mw 7.0 event before 1952, which is consistent with the actual occurrence of the M 6.8 earthquake in 1952 (Fig. 3B). After accounting for the stress perturbation of the 1952 earthquake (“Methods” section), we assess the $G_0$, $G_c$ and the maximum possible earthquake magnitude in the current stress state (Fig. 3C, D). Because the 1952 earthquake has released part of the stored elastic energy, the asperity AZF-3 has less possibility to produce a Mw 7 event in its current state. But its neighboring asperities can still

accommodate Mw > 7 events and their maximum magnitudes slightly increase due to the loading from the 1952 rupture. Due to the lack of direct constraints on the 1952 M 6.8 earthquake in the current model, the possibility of a future Mw 7.5 event bridging AFZ-3 and AFZ-4, similar to the 1536 M 7.5 event, cannot be excluded and thus is worthy of further investigation. On the AZF, we primarily focus on the Anninghe fault and the northern segment of the Zemuhe fault, as the shear stress on the southern segment of the Zemuhe fault is currently insufficient to accommodate large events. We also use this theory to estimate the maximum possible earthquake magnitude on the DLSF, although its coupling is less finely resolved than that of the AZF. Due to the absence of paleoseismic evidence for the DLSF, we assume that the elapsed time since the last large earthquake is either 500 or 1000 years, to obtain a lower bound estimate of the maximum possible earthquake magnitude. We find the maximum possible earthquake magnitude on the DLSF could also reach Mw 7 (Fig. S3) if the elapsed time since the last large earthquake is longer than 500 years.

To confirm these theoretical predictions, we first set a 3D dynamic rupture model with a vertical planar fault intersecting the free surface in a linear elastic medium. In the numerical model, we implement the power and exponential slip-weakening friction laws, respectively, whose coefficients are determined based on the data of global earthquakes within the relevant range of slip values (“Methods” section). These two friction laws represent two end-member regimes of thermal pressurization³⁸ caused by the shear heating of fault zone pore fluids. The exponential friction law is also proposed in laboratory experiments³⁹. The current shear stress on the fault has been building up since the latest large earthquake, during which the frictional strength dropped to a residual value $\tau_r$. Due to the lack of direct constraints, we need to indirectly infer $\tau_r$ value of the last big event. For the exponential friction, we assume $\tau_r$ equals the minimum dynamic strength $\tau_d$ (“Methods” section). Note that the absolute value of $\tau_d$ is irrelevant in elastic rupture dynamics on vertical planar faults⁹. For the power-law friction, we assume $\tau_r$ is larger than $\tau_d$, corresponding to an empirical slip of the last big event (“Methods” section). As supported by laboratory experiments of large slips³⁹, the assumption of $\tau_r = \tau_d$ for the exponential friction relies on the estimated slip of the last big event being sufficiently large. The assumption for the power friction relies on the estimated slip of the historical big event: overestimated slips lead to underestimated $\tau_r$ and thus result in
<figure><img src="image_7.png" /><figcaption>**Fig. 4** | **Effects of fault geometries on maximum possible magnitude of earthquakes on the Anninghe-Zemuhe fault (AZF).** A The surface traces of curved and planar faults (legend) used in dynamic rupture simulations. B The maximum possible earthquake magnitude in 2023 as a function of hypocenter location predicted by simulations with curved and planar faults (legend in A). The simulation results</figcaption></figure>

are compared with the theoretical prediction. C Accumulated shear stress on curved fault (top) and superimposition of slip for four earthquake scenarios on curved (upper middle), smoothed (lower middle), and planar (bottom) faults, each assuming a different nucleation location (marked by purple stars).

underestimated potential earthquake sizes. Despite these different assumptions of shear stresses and friction laws, we find that the theoretical and numerical predictions (Figs. 3 and S12) are quantitatively consistent in different rupture scenarios within the uncertainties in fitting the $G_c$ scaling data. Both the theory and numerical model predict that the maximum possible earthquake magnitude may increase to Mw > 7.3 on the AZF if the earthquake occurred in 2200 (Figs. 3D and S13).

We then conduct additional 3D dynamic rupture models with a curved fault as used in the coupling inversion model (red curve in Fig. 4A), characterized by a significant bending within segment AZF-4. For comparison, we also consider a model with smoothed fault trace (blue curve in Fig. 4A). Based on the current stress state of the AZF constrained in this study, we find that fault bending is not the primary factor controlling earthquake sizes on the AZF, and that the numerical simulations are quantitatively consistent with the theoretical prediction regardless of different fault geometries on the AZF (Fig. 4B, C). The explanation may be that the dominant wavelength of fault curvature on the AZF is larger than the length scales of rupture depths or seismic asperities. Consequently, the first-order rupture dynamics on the AZF can be well predicted by 3D dynamic fracture mechanics based on a planar fault. While in the future, it is necessary to consider the effects of fault complexities, such as step-overs on the DLSF, despite its observations being less constrained compared to those of the AZF.

## Discussion

In earthquake cycle models, small seismicity can be seen as failed initiations of potential large earthquakes⁴⁰. The seismicity in asperity AZF-1 is more active than in the other asperities (Fig. 3), which implies that this asperity has more frequent chances to initiate a large event with maximum possible magnitude reaching up to Mw 7. To the north of this segment, the Mw 6.6 Luding earthquake occurred in 2022 (ref. 41), triggered many aftershocks⁴², and loaded this fault segment by 0.2 MPa (Fig. S4). Based on these seismic observations and our models, we suggest that asperity AZF-1 has higher seismic hazard than the other segments and is currently ready to host a Mw 7 event.

Due to the high-coupling ratio and long elapsed times, the current accumulated slip deficits on asperities AZF-1 and AZF-4 reach ~3 m, close to the coseismic slip of historical events as revealed by paleoseismic investigations⁴³. The slip deficit on asperities AZF-5 and AZF-6

is less than 1.0 m due to the short time elapsed since their last ruptures in 1732 and 1850. The accumulated slip deficit is less than 0.5 m on asperity AZF-3 due to its recent last rupture in 1952. We estimate the recurrence intervals of big events on these asperities based on 3D dynamic rupture theory ("Methods" section): $T = s_b/S$, where S and s_b are the critical slip for recurred ruptures and the slip deficit rate derived from the kinematic inversion, respectively. We estimate S at any given fault point as the slip deficit value for which the point reaches the necessary condition for runaway rupture, $G_0 \approx G_c$. Because $G_0$ increases faster than $G_c$ with the accumulated slip deficit, S can be uniquely determined ("Methods" section). Note that this estimate only gives a lower bound on the rupture time because the earthquake may nucleate later when the runaway condition is satisfied. As shown in Fig. 2D, the estimated return periods are about 500 years on these asperities, quantitatively consistent with that inferred from geological studies⁴³,⁴⁴. Two sampling sites located in the low coupling zones (Ganhaizi and Daqingliangzi) have long return periods, which is also close to our estimate on the low coupling patches.

In summary, we combine geodetic, seismological, and geological observations with physics-based models to assess the seismic hazard of the Anninghe-Zemuhe fault system. While faulting is also affected by additional complicated processes whose parameters are poorly constrained currently, such as fault zone fluids, geometric complexities, and the state of absolute stress, we demonstrated that a useful estimate of the maximum magnitude of potential earthquakes is readily available based on the accumulated slip deficit on faults and dynamic rupture models. Improved constraints on fault coupling and stress state, thanks to future developments of geodetic and seismological observations, together with advances in dynamic fracture theory and numerical simulations like the ones demonstrated here, should enable a more reliable, time-dependent evaluation of the location and magnitude of potential large earthquakes.

## Methods

### Inversion of fault coupling ratio

The fault coupling models of the AZF, as shown in previous studies²⁶⁻²⁸, may suffer from the sparsely distributed GNSS stations near the fault (Fig. 1). To better constrain the coupling ratio distribution along the AZF, we built 24 new near-fault GNSS stations in 2018 and surveyed them three times in 2019, 2020 and 2022 (Fig. 1A). In order to reduce
seasonal effects, field surveys were performed at the similar time of the year (July). The occupation of each survey was larger than 72 h to ensure the stability of the observation. We solve the velocities on these new GNSS stations using the strategy shown in Diao et al.⁴⁵ and combine them with the existing GNSS velocities⁴⁶ for the kinematic inversion following the method of Wang and Shen⁴⁶. The solved GNSS velocities and associated uncertainties on these new stations are shown in Table S1. To model the crustal deformation and invert for the fault coupling ratio, a combined 3D deformation model as shown in Eq. (1) has been widely used²⁶,

$$V_{\text{obs}} = V_B + V_{\varepsilon} + V_p + V_{\text{backslip}} \quad (1)$$

where $V_{\text{obs}}$ are the observed GNSS velocities, $V_B$ and $V_{\varepsilon}$ the velocities due to steady state block rotation and interior strain²⁶, respectively, $V_p$ the postseismic viscoelastic relaxation effect of large historical earthquakes, and $V_{\text{backslip}}$ the deformation caused by fault coupling.

$V_p$ is accounted for in our model by simulating the postseismic viscoelastic effect of historical earthquakes since 1480. For this purpose, we use the PSGRN/PSCMP code⁴⁷ with a layered infinite half-space model that consists of an elastic upper crust and two underlying viscoelastic layers representing the lower crust and upper mantle (Fig. S5). The elastic and viscoelastic parameters in the model are from ref. 48 and refs. 45,49, respectively. The postseismic viscoelastic relaxation effect is driven by the coseismic rupture of historical earthquakes, of which the parameters are inferred from refs. 24,50–52 and listed in Table S2. As shown in Fig. S6, this effect remains less than 0.5 mm/a despite different viscosities being used. Beside the viscoelastic relaxation effect due to stress loading of historical events ($V_p$ in Eq. (1)), we also consider the viscoelastic effect caused by the interseismic stress accumulation due to the fault coupling. Based on the theoretical derivation shown in ref. 53, we consider the interseismic viscoelastic effect by constructing the viscoelastic Green’s function of $G(x,t) = [E(x) + \text{Vis}(x,t)H(t)]$ for slip deficit inversions (Eq. (2)), where $H(t)$ is the Heaviside function, $E(x)$ the elastic displacement, and $\text{Vis}(x,t)$ the time-dependent displacement induced by the interseismic viscoelastic relaxation. Considering the long earthquake recurrence period on the AZF, we calculate $G(x,t)$ by simulating a fully relaxed viscoelastic effect ($t \to \infty$ using the layered viscoelastic model (Fig. S5).

$V_B$ and $V_{\varepsilon}$ are calculated by using an Euler vector and a strain rate tensor⁵⁴ that are solved by the least square method from the GNSS velocities within each block. We divide the study area into three blocks (Model 1 in Fig. S7) based on the strategy inferred from an Euler pole clustering technique⁵⁵. We also test the other block division strategy (Model 2 in Fig. S7) to investigate whether the LJ-XJHF shall be defined as the block boundary in the model. In Model 1 and Model 2 we use different coverages of GNSS stations to define the motion of blocks west of the AZF, i.e. Chuandian block and Dianzhong block in Fig. S7. Modeling results indicate that Model 2 leads to a slip rate of 9.2–5.7 mm/a along the AZF (Model 2 in Fig. S7) that is higher than the geological result of 6.5 ± 1.0 mm/a³¹. On the contrary, Model 1 that neglects the XJHF, generates a slip rate of AZF (6.1–4.9 mm/a) that agrees better with the geologic investigations.

After decomposing $V_B$, $V_{\varepsilon}$, and $V_p$ from the observations, the residual deformation is assumed to be caused by the fault coupling that is represented by $V_{\text{backslip}}$. We set the fault trace based on the geologic data⁵⁶ and assume the AZF to be vertical following previous studies²⁶⁻²⁸,⁵⁷. We further subdivide the fault into discrete patches with a size of 5 × 5 km, and invert for back-slip rate utilizing a constrained least squares method⁴⁵. To make the solution stable, we add an a-priori smoothing constraint in the cost function,

$$F(s_b) = \|Gs_b - y\|^2 + \alpha^2 \|Hs_b\|^2 \quad (2)$$

where $s_b$ is the back-slip rate vector on fault patches, G the aforementioned viscoelastic Green’s function, y the deformation component related to fault back-slip rate (in Eq. (1)), H the finite difference approximation of the Laplacian operator, and α the smoothing factor that controls the trade-off between model roughness and data misfit (Fig. S8). After obtaining the back-slip rate, we obtain the fault coupling ratio (r) by $r = s_b/s$, where s is the fault slip rate. We estimate the uncertainty of the inverted slip deficit rate and coupling ratio using the statistic method⁵³. The details are as follows: we first generate a set of Gaussian noise signals using the observation uncertainty as the standard deviation. Then, we add them to the observation data to get several sets of disturbed synthetic data. For each of the disturbed datasets, we repeat the inversion process to get the coupling ratio on each fault patch and the mean coupling ratio above 15 km. These steps were run 500 times in order to estimate the uncertainty of the parameters statistically.

Moreover, to test the effect of secondary faults within the Chuandian block (Lijiang-Xiaojinhe fault and Litang fault), we forward calculate the crustal deformation due to the motion of these secondary faults using the “buried-dislocation” model⁵⁴ and associated fault slip rate⁵⁸,⁵⁹. As shown in Fig. S9, such effect is ~0.5 mm/a near the AZF. After removing this effect from the GNSS velocities, we re-run Model 1 to get the fault coupling distribution (Fig. S9). Not surprisingly, we find that the coupling model remains stable and the effect of these secondary faults is negligible.

## Estimates of the b value

Based on the seismic records in the past 14 years (2009–2023), we update the earthquake catalog of the study area using the TomoDD code⁶⁰, and obtain the precision location of more than 94400 ML > 1.5 events in this region. Based on this refined earthquake catalog (Supplementary Data 1), we estimate the b values homogeneously over space by employing the maximum-likelihood method within the ZMAP software⁶¹, with a dense spatial grid (0.03 × 0.03°) and a sampling volume of circular shape ($r = 25$ km). We also estimate the magnitude of completeness (Mc) on each grid using the Goodness-of-Fit test⁶². For a good fit to the G-R law⁶¹, we only calculate the b values on the grids where the number of events (magnitude > Mc) exceeds $N_{\text{min}}$. We collect the b values on grids near the AZF (<5 km) to build the profile shown in Fig. 2C. To verify the robustness of the inferred b values, we conduct additional tests. First, we increase $N_{\text{min}}$ from 30 to 60 and observe the resulting variation in the calculated b values (Fig. S10). Then, we further increase $N_{\text{min}}$ to 80 and 100, but with an enlarged sampling radius $r$ of 40 km (Fig. S11). We find that the primary features of the inferred b values are qualitatively consistent, despite low $N_{\text{min}}$ values resulting in large uncertainties and high $N_{\text{min}}$ values leading to a highly smoothed variation (Figs. S10, S11). In this study, we therefore choose $N_{\text{min}} = 50$ with a sampling radius $r$ of 25 km, following ref. 37, to ensure sufficient spatial sampling while preserving spatial variations of the b values in each asperity. In addition, the frequency-magnitude distributions within selected patches also show that the estimates of the b value are robust (Fig. S10).

## Scaling relations between fracture energy and slip

To constrain frictional properties on faults, we compile the datasets of global earthquakes from the literature⁶³⁻⁶⁶ and determine the scaling relations between fracture energy and final slip (Fig. S12A). It has been proposed that fracture energy $G_c$ increases with increasing final slip D as inferred in natural faults³⁸,⁶⁷,⁶⁸ and also observed in laboratory experiments⁶⁹. As a first-order approximation derived from off-fault inelastic dissipation⁷⁰,⁷¹ and thermal pressurization³⁸, the scaling relation between the fracture energy and final slip can be expressed as a power-law, $G_c \approx BD^n$, where B and n are coefficients to be determined. Alternatively, laboratory experiments³⁹ show that the friction weakening behaviors of natural rocks can be described by an exponential
slip-weakening friction, $\tau = (\tau_s - \tau_d)e^{-D/d_c} + \tau_d$, where $(\tau_s - \tau_d)$ and $d_c$ are the strength drop and the critical slip-weakening distance. For this friction, the scaling relation between fracture energy and final slip is $G_c = (\tau_s - \tau_d)d_c \cdot (1 - (1+D/d_c)e^{-D/d_c})$, where $(\tau_s - \tau_d)$ and $d_c$ are coefficients to be determined. In this study, we separately constrain these two sets of free coefficients and use them for the seismic hazard assessment (Fig. 3). The best fitting coefficients are $B=1.8$ and $n=0.85$ for the power-law slip-weakening friction, and $\tau_s - \tau_d = 8$ MPa and $d_c = 0.2$ m for the exponential slip-weakening friction (Fig. S12A). The two friction laws are comparably consistent with the fracture energy scaling data within the range of slip values that are most relevant for this study (0.1–1 m). To do the sensitivity analysis for the theoretical models, we vary the value of $B$ by a factor of 1.4 (Fig. 3). Note that in this study we only include the data of $G_c$ constrained by dynamic models in the literature$^{63–66}$ for our fitting of friction coefficients, therefore, our constrained value of $B$ is different from that in ref. 9, which also included the experimental data in their fitting.

## Dynamic rupture simulations

We set 3D dynamic rupture simulations on a vertical planar or curved fault of 20 km width intersecting the free surface in a linear elastic and homogeneous medium. The computational domain has a size of 200 × 100 × 70 km, which is large enough to avoid contamination by waves reflected from the artificial boundaries during the rupture times. The P wave speed, S wave speed, and density of the medium of the elastic model are uniformly assumed as 5770 m/s, 3330 m/s, and 2705 g/m$^3$ based on the parameters in the upper layers shown by Liu et al.$^{48}$. As the variation of shear modulus is less than 10% within the seismogenic zone$^{48}$, the specific choice of the velocity models has a secondary effect on the numerical simulations. We use the spectral element software SPECFEM3D$^{72}$ for the dynamic simulation. The time step is set to 0.007 s (Table S3), meeting the Courant-Friedrichs-Lewy stability condition, a criterion widely used to ensure the stability of the numerical simulation. We set a grid size of 250 m to guarantee sufficient numerical resolution.

We use both the power and exponential slip-weakening friction laws in the numerical simulations. To confirm the theoretical model, we first use the nonlinear friction model proposed by ref. 73: $\tau = \tau_d + (\tau_s - \tau_d)/(1 + \frac{\delta}{pd_c})^p$, where $p$ is the power-law coefficient, $\tau_s$, $\tau_d$, and $d_c$ are the static strength, the dynamic strength, and the critical slip-weakening distance, respectively. In the simulations, we set $p=1-n=0.15$, $d_c=0.01m$, and $\tau_s - \tau_d = 28MPa$, as determined by the data of global earthquakes within the relevant range of slip values (red curve in Fig. S12A). This friction law produces a similar scaling relation as that used in the theoretical model (Fig. S12A). The initial shear stress is set as $\tau_{ini} = \Delta\tau_{accu} + \tau_r$, where we assume that the latest large earthquakes dropped the shear stresses to the minimum strength that is corresponding to a slip of 5 m, thus we set $\tau_r = \tau_d + (\tau_s - \tau_d)/(1 + \frac{\delta}{pd_c})^p$. For the exponential slip-weakening friction law, we assume $\tau = (\tau_s - \tau_d)e^{-\delta/d_c} + \tau_d$, where $\tau_s$ and $\tau_d$, and $d_c$ are the static strength, the dynamic strength, and the critical slip-weakening distance, respectively. The friction coefficients are also determined based on the data of observed global earthquakes by fitting the scaling of the fracture energy (blue curve in Fig. S12A), where $\tau_s - \tau_d = 8$ MPa and $d_c = 0.2m$. Note that the absolute value of $\tau_d$ is not important in these two friction laws on vertical planer faults$^{70}$, therefore here we artificially set $\tau_d = 1.2$ MPa. In the exponential friction law, the initial shear stress is set as $\tau_{ini} = \Delta\tau_{accu} + \tau_d$, where we assume that the latest large earthquakes dropped the shear stresses to the minimum dynamic strength $\tau_d$ and the current shear stresses built up since then. We list the key parameters used for dynamic rupture simulations in Table S3.

To estimate $\Delta\tau_{accu}$, we first use the coupling ratio to calculate the shear stress changes per year $\tau$ and then calculate the total

accumulated shear stress based on the elapsed time $T$ since the latest large earthquakes, $\Delta\tau_{accu} = \tau T$ (Fig. 2B). For the simulations with stresses at the current time (Fig. 3D), we account for the stress perturbation caused by the 1952 earthquake. As there is no slip model available for that event, we use the stress change from a numerical model that produces a similar magnitude earthquake.

We use a time-dependent weakening to nucleate the ruptures, where the nucleation zone is 6 km and the nucleation speed is 1200 m/s. Rupture propagation is spontaneously controlled by the slip-weakening friction law outside the nucleation zone. In the simulations, the depth of the nucleation is 6 km and the horizontal nucleation location is varied along strike in different models to search for the maximum possible earthquake magnitude as a function of the nucleation location. For all the ruptures, we calculate the final seismic moment and compare them with the theoretical prediction (Fig. 3).

## Theoretical model based on 3D dynamic fracture mechanics

The rupture arrest locations and earthquake sizes are estimated based on the theory of 3D dynamic fracture mechanics on long faults with finite rupture width$^{9,10}$. The 3D theory$^{9,10}$ suggests a simple criterion for determining the arrest location of dynamic ruptures by evaluating the rupture potential: $\phi = \int_{x_0}^{x_0+L} (1 - G_c(x)/G_0(x))dx = 0$, where $L$ is predicted along-strike arrest location distance from the hypocenter (located at $x_0$), $G_0(x)$ and $G_c(x)$ are the potential elastic energy and dissipated fracture energy as a function of the along-strike location $x$, respectively. The criterion applies to each of the rupture fronts on bilateral ruptures. The potential elastic energy $G_0(x)$ is calculated based on the seismic coupling map and the elapsed time since the last M-7.5 large earthquakes. To obtain $G_0(x)$, we first estimate the total slip deficit $D_{accu}$ based on the elapsed time $T$ since the latest large earthquakes, and then calculate the accumulated shear stress $\Delta\tau_{accu}$ on fault. The density of the potential elastic energy accumulated on the seismogenic fault is defined by $0.5\Delta\tau_{accu}D_{accu}$. Following previous method$^{10}$, we average $0.5\Delta\tau_{accu}D_{accu}$ across the seismogenic depth (i.e., $W = 20km$ in this paper) to obtain $G_0(x)$ along strike (Fig. 3A, C). The fracture energy $G_c$ is estimated based on potential earthquake slip $D$ and the scaling relation $G_c(x) = BD^n$ with $B=1.8$ and $n=0.85$ from the best fitting of the observations. For the models in the current time (Fig. 3D), we account for the 1952 earthquake by subtracting the approximated slip of this earthquake from the total slip deficit. We then calculate the density of the potential fracture energy by using $BD_{accu}^n$ and average it across the seismogenic depth to obtain $G_c(x)$ along strike. Finally, we use the criterion of rupture potential to determine the final rupture length $L$ given an assumed nucleation location and compute the final seismic moment by integrating the released slip along the ruptured area.

## Data availability

The velocity of newly-developed GNSS stations is listed in Table S1, and the other GNSS velocities can be found in the supplementary files of ref. 46. The refined seismic catalog used for $b$ value calculation is shown in Supplementary Data 1.

## Code availability

The PSGRN/PSCMP code for postseismic deformation modeling and viscoelastic Green's function calculation, as well as the SDM code for the slip deficit inversions, are available at the public server of German Research Centre for Geosciences (GFZ): ftp://ftp.gfz-potsdam.de/pub/home/turk/wang/. The open-source software SPECFEM3D used in dynamic rupture simulations is available at Zenodo (https://doi.org/10.5281/zenodo.10413988). The open-source software ZMAP used for $b$ value calculations is available at: https://doi.org/10.5281/zenodo.3470871 or https://github.com/zmap/zmap. We drew the figures using the open-source software GMT5 (http://gmt.soest.hawaii.edu/projects/gmt).
## References

1. Brace, W. F. et al. Stick-slip as a mechanism for earthquakes. *Science* **153**, 990–992 (1966).

2. Bakun, W. H. & Lindh, A. G. The Parkfield, California, earthquake prediction experiment. *Science* **229**, 619–624 (1985).

3. Chen, T. & Lapusta, N. On behaviour and scaling of small repeating earthquakes in rate and state fault models. *Geophys. J. Int.* **218**, 2001–2018 (2019).

4. Nocquet, J.-M. et al. Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake. *Nat. Geosci.* **10**, 145–149 (2017).

5. Jackson, D. D. & Kagan, Y. Y. The 2004 Parkfield earthquake, the 1985 prediction, and characteristic earthquakes: lessons for the future. *Bull. Seismol. Soc. Am.* **96**, S397–S409 (2006).

6. Barbot, S. et al. Under the hood of the earthquake machine: toward predictive modeling of the seismic cycle. *Science* **336**, 707–710 (2012).

7. Kaneko, Y. & Lapusta, N. Variability of earthquake nucleation in continuum models of rate-and-state faults and implications for aftershock rates. *J. Geophys. Res. Solid Earth* **113**, B12312 (2008).

8. Yang, H. et al. Deriving rupture scenarios from interseismic locking distributions along the subduction megathrust. *J. Geophys. Res. Solid Earth* **124**, 10376–10392 (2019).

9. Weng, H. & Ampuero, J. P. The dynamics of elongated earthquake ruptures. *J. Geophys. Res. Solid Earth* **124**, 8584–8610 (2019).

10. Weng, H. & Ampuero, J.-P. Continuum of earthquake rupture speeds enabled by oblique slip. *Nat. Geosci.* **13**, 817–821 (2020).

11. Metois, M. et al. Interseismic coupling, segmentation and mechanical behavior of the central Chile subduction zone. *J. Geophys. Res. Solid Earth* **117**, B03406 (2012).

12. Jolivet, R. et al. Aseismic slip and seismogenic coupling along the central San Andreas Fault. *Geophys. Res. Lett.* **42**, 297–306 (2015).

13. Stevens, V. L. & Avouac, J. P. Interseismic coupling on the main Himalayan thrust. *Geophys. Res. Lett.* **42**, 5828–5837 (2015).

14. Hok, S. et al. Dynamic rupture scenarios of anticipated Nankai-Tonankai earthquakes, southwest Japan. *J. Geophys. Res.* **116**, B12319 (2011).

15. Harris, R. A. et al. Fault steps and the dynamic rupture process: 2-D numerical simulations of a spontaneously propagating shear fracture. *Geophys. Res. Lett.* **18**, 893–896 (1991).

16. Harris, R. A. & Day, S. M. Dynamics of fault interaction: parallel strike-slip faults. *J. Geophys. Res.* **98**, 4461–4472 (1993).

17. Ulrich, T. et al. Dynamic viability of the 2016 Mw 7.8 Kaikoura earthquake cascade on weak crustal faults. *Nat. Commun.* **10**, 1213 (2019).

18. Ramos, M. D. et al. Assessing margin-wide rupture behaviors along the Cascadia megathrust with 3-D dynamic rupture simulations. *J. Geophys. Res.* **126**, e2021JB022005 (2021).

19. Duan, B. & Oglesby, D. D. Multicycle dynamics of nonplanar strike-slip faults. *J. Geophys. Res.* **110**, B03304 (2005).

20. Dieterich, J. H. & Richards-Dinger, K. B. Earthquake recurrence in simulated fault systems. *Pure Appl. Geophys.* **167**, 1087–1104 (2010).

21. Wesnousky, S. G. Predicting the endpoints of earthquake ruptures. *Nature* **444**, 358–360 (2006).

22. Biasi, G. P. & Wesnousky, S. G. Bends and ends of surface ruptures. *Bull. Seismol. Soc. Am.* **107**, 2543–2560 (2017).

23. Shaw, B. E. et al. A physics-based earthquake simulator replicates seismic hazard statistics across California. *Sci. Adv.* **4**, eaau0688 (2018).

24. Wen, X. Z. et al. Historical pattern and behavior of earthquake ruptures along the eastern boundary of the Sichuan-Yunnan faulted-block, southwestern China. *Phys. Earth Planet. Inter.* **168**, 16–36 (2008).

25. Sun, H. et al. Paleoearthquake history along the southern segment of the Daliangshan fault zone in the southeastern Tibetan Plateau. *Tectonics* **38**, 2208–2231 (2019).

26. Jiang, G. et al. Geodetic imaging of potential seismogenic asperities on the Xianshuihe-Anninghe-Zemuhe fault system, southwest China, with a new 3-D viscoelastic interseismic coupling model. *J. Geophys. Res. Solid Earth* **120**, 1855–1873 (2015).

27. Li, Y. et al. Heterogeneous interseismic coupling along the Xian-shuihe-Xiaojiang fault system, eastern Tibet. *J. Geophys. Res. Solid Earth* **126**, e2020JB021187 (2021).

28. Li, Y. et al. Interseismic coupling, asperity distribution, and earthquake potential on major faults in southeastern Tibet. *Geophys. Res. Lett.* **50**, e2022GL101209 (2023).

29. Wang, K. et al. Deformation cycles of subduction earthquakes in a viscoelastic Earth. *Nature* **484**, 327–332 (2012).

30. Li, S. et al. Revisiting viscoelastic effects on interseismic deformation and locking degree: a case study of the Peru-North Chile subduction zone. *J. Geophys. Res. Solid Earth* **120**, 4522–4538 (2015).

31. Xu, X. W. et al. Newly tectonic deformation and its dynamic source of the Sichuan-Yunnan Block. *Sci. China Earth Sci.* **46**, 210–226 (2003).

32. He, H. et al. Newly-generated Daliangshan fault zone—Shortcutting on the central section of Xianshuihe-Xiaojiang fault system. *Sci. China Earth Sci.* **51**, 1248–1258 (2008).

33. Bletery, Q. et al. Mega-earthquakes rupture flat megathrusts. *Science* **354**, 1027–1031 (2016).

34. Manighetti, I. et al. Earthquake scaling, fault segmentation, and structural maturity. *Earth Planet. Sci. Lett.* **253**, 429–438 (2007).

35. Xu, S. et al. Fault strength and rupture process controlled by fault surface topography. *Nat. Geosci.* **16**, 94–100 (2023).

36. Schorlemmer, D., Wiemer, S. & Wyss, M. Variations in earthquake-size distribution across different stress regimes. *Nature* **437**, 539–542 (2005).

37. Nanjo, K. Z. & Yoshida, A. A b map implying the first eastern rupture of the Nankai Trough earthquakes. *Nat. Commun.* **9**, 1117 (2018).

38. Viesca, R. C. & Garagash, D. I. Ubiquitous weakening of faults due to thermal pressurization. *Nat. Geosci.* **8**, 875–879 (2015).

39. Di Toro, G. et al. Fault lubrication during earthquakes. *Nature* **471**, 494–498 (2011).

40. Lapusta, N. & Rice, J. R. Nucleation and early seismic propagation of small and large events in a crustal earthquake model. *J. Geophys. Res. Solid Earth* **108(B4)**, 2205 (2003).

41. Li, Y. et al. Coseismic slip model of the 2022 Mw 6.7 Luding (Tibet) earthquake: pre- and post-earthquake interactions with surrounding major faults. *Geophys. Res. Lett.* **49**, e2022GL102043 (2022).

42. Zhang, L. et al. 2022 Mw 6.6 Luding, China, earthquake: a strong continental event illuminating the Moxi seismic gap. *Seismol. Res. Lett.* **94**, 2129–2142 (2023).

43. Wang, H. et al. A 3400-year-long paleoseismologic record of earthquakes on the southern segment of Anninghe fault on the southeastern margin of the Tibetan Plateau. *Tectonophysics* **628**, 206–217 (2014).

44. Wang, H. et al. Paleoearthquakes on the Anninghe and Zemuhe fault along the southeastern margin of the Tibetan Plateau and implications for fault rupture behavior at fault bends on strike-slip faults. *Tectonophysics* **721**, 167–178 (2017).

45. Diao, F. et al. Fault behavior and lower crustal rheology inferred from the first seven years of postseismic GPS data after the 2008 Wenchuan earthquake. *Earth Planet. Sci. Lett.* **495**, 202–212 (2018).

46. Wang, M. & Shen, Z. Present-day crustal deformation of continental China derived from GPS and its tectonic implications. *J. Geophys. Res. Solid Earth* **125**, e2019JB018774 (2020).

47. Wang, R. et al. PSGRN/PSCMP—a new code for calculating co- and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory. *Comput. Geosci.* **32**, 527–541 (2006).

48. Liu, Y. et al. The high-resolution community velocity model V2.0 of southwest China, constructed by joint body and surface wave
tomography of data recorded at temporary dense arrays. *Sci. China Earth Sci.* **66**, 2368-2385 (2023).

49. Huang, M. H. et al. Probing the lithospheric rheology across the eastern margin of the Tibetan Plateau. *Earth Planet. Sci. Lett.* **396**, 88-96 (2014).

50. Papadimitriou, E. et al. Earthquake triggering along the Xianshuihe fault zone of Western Sichuan, China. *Pure Appl. Geophys.* **161**, 1683-1707 (2004).

51. Wen, X. Z. et al. Tectonic dynamics and correlation of major earthquake sequences of the Xiaojiang and Qujiang-Shiping fault systems, Yunnan, China. *Sci. China Earth Sci.* **54**, 1563-1575 (2011).

52. Zhu, H. & Wen, X. Z. Effects of static stress triggering of the major earthquake sequence on Xiaojiang-Zemuhe fault zone. *Earth Sci.* **37**, 199-206 (2012).

53. Diao, F. et al. Slip Rate Variation Along the Kunlun Fault (Tibet): results from new GPS observations and a viscoelastic earthquake-cycle deformation model. *Geophys. Res. Lett.* **46**, 2524-2533 (2019).

54. Savage, J. C. et al. Strain accumulation and rotation in the eastern California shear zone. *J. Geophys. Res.* **106**, 21995-22007 (2001).

55. Xu, R. & Stamps, D. S. Strain accommodation in the Daliangshan Mountain area, southeastern margin of the Tibetan Plateau. *J. Geophys. Res. Solid Earth* **124**, 9816-9832 (2019).

56. Deng, Q. et al. Basic characteristics of active tectonics of China. *Sci. China Earth Sci.* **46**, 356-372 (2003).

57. Yao, S. & Yang, H. Hypocentral dependent shallow slip distribution and rupture extents along a strike-slip fault. *Earth Planet. Sci. Lett.* **578**, 117296 (2022).

58. Xu, X. et al. Average slip rate, earthquake rupturing segmentation and recurrence behavior on the Litang fault zone, western Sichuan Province, China. *Sci. China Earth Sci.* **48**, 1183-1196 (2005).

59. Chevalier, M. L. et al. Tectonic-geomorphology of the Litang fault system, SE Tibetan Plateau, and implication for regional seismic hazard. *Tectonophysics* **682**, 278-292 (2016).

60. Zhang, H. & Thurber, C. Development and applications of double-difference seismic tomography. *Pure Appl. Geophys.* **163**, 373-403 (2006).

61. Wiemer, S. A software package to analyze seismicity: ZMAP. *Seismol. Res. Lett.* **72**, 373-382 (2001).

62. Gutenberg, B. & Richter, C. F. Frequency of earthquakes in California. *Bull. Seismol. Soc. Am.* **34**, 185-188 (1944).

63. Mai, P. et al. On scaling of fracture energy and stress drop in dynamic rupture models: consequences for near-source ground-motions. *Earthquakes Radiat. Energy Phys. Faulting* **170**, 283-293 (2006).

64. Gallovič, F. et al. Bayesian dynamic finite-fault inversion: 2. Application to the 2016 Mw 6.2 Amatrice, Italy, earthquake. *J. Geophys. Res. Solid Earth* **124**, 6970-6988 (2019).

65. Weng, H. & Yang, H. Constraining frictional properties on fault by dynamic rupture simulations and near-field observations. *J. Geophys. Res. Solid Earth* **123**, 6658-6670 (2018).

66. Yao, S. & Yang, H. Rupture dynamics of the 2012 Nicoya Mw7.6 earthquake: evidence for low strength on the megathrust. *Geophys. Res. Lett.* **47**, e2020GL087508 (2020).

67. Rice, J. Heating and weakening of faults during earthquake slip. *J. Geophys. Res. Solid Earth* **111**, B05311 (2006).

68. Tinti, E. et al. Estimates of slip weakening distance for different dynamic rupture models. *Geophys. Res. Lett.* **31**, L02611 (2004).

69. Nielsen, S. et al. G: fracture energy, friction and dissipation in earthquakes. *J. Seismol.* **20**, 1187-1205 (2016).

70. Andrews, D. J. Rupture dynamics with energy loss outside the slip zone. *J. Geophys. Res. Solid Earth* **110**, B01307 (2005).

71. Gabriel, A. A. et al. Source properties of dynamic rupture pulses with off-fault plasticity. *J. Geophys. Res. Solid Earth* **118**, 4117-4126 (2013).

72. Ampuero, J. P. et al. Nucleation of rupture under slip dependent friction law: simple models of fault zone. *J. Geophys. Res. Solid Earth* **107**, ESE-2 (2002).

73. Chambon, G. et al. Frictional response of a thick gouge sample: 2. Friction law and implications for faults. *J. Geophys. Res. Solid Earth* **111**, B09309 (2006).

## Acknowledgements

This work was supported by the National Key Research and Development Program of China (no.2023YFC3007302, F.D., H.W., and Z.S.) and the National Natural Science Foundation of China (no. 42374066, H.W., no. 41731072, X.X. and F.D.). We thank the survey crews for GNSS data collection in the fieldwork, and Yuebing Wang for assistance in GNSS data processing.

## Author contributions

F.D. and H.W. designed the study. F.D. carried out the GNSS data collection and processing, fault coupling inversion, and results analysis. H.W. performed the numerical experiments and theoretical analysis. F.L. carried out the earthquake catalog refinement and *b* value calculation. F.D. and H.W. wrote the initial manuscript. J.P.A., Z.S., R.W., and X.X. discussed the results and contributed to writing the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-024-51313-w.

**Correspondence** and requests for materials should be addressed to Faqi Diao or Huihui Weng.

**Peer review information** *Nature Communications* thanks the anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2024