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

Domain: earth_science
Published: 2024-08-12
Editorial outcome: minor_revision

Abstract
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 Mw 7.0 earthquakes with maximum magnitudes increasing to Mw 7.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.

Introduction
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 fault1,2 and small repeating earthquakes3. 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 supercycle4. Predicting the recurrence times of large earthquakes is challenging, even for the simple characteristic model5,6, mainly due to the unpredictable nucleation conditions7 and the uncertainties in stress states8. While we still cannot predict where and when earthquakes will start, recent advances9,10 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 models11–13 or dynamically-based models8,14–20 on various fault systems, including estimates of their slip budget, have been used to infer their potential seismic moments. Specifically, fault geometric complexities21,22, such as fault bends and step-overs, have been integrated into single rupture simulations15–18 and earthquake cycle models19,20. Furthermore, complex fault systems have also been incorporated with a physics-based statistical simulator23 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 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 events24, whereas there is barely paleoseismic evidence of large events on the DLSF in the past 500 years25. 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 earthquakes26–28. 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 limit 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 areas29. 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 ratio30. 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/a31. 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/a32. 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 maturities32: the mature AZF has more continuous surface traces than the immature DLSF (Fig. 1A). Seismological33,34 and experimental34 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 addition, lacking paleoseismic evidence of large earthquakes on the immature DLSF25 is qualitatively consistent with previous studies of fault roughness33,35. 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 comparison 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 accumulated 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 historical earthquakes in the past ~500 years, delineated by comprehensive analysis of multi-disciplinary observations24. 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 level36,37. 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, 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 theory9, rupture propagation and arrest on a fault are controlled by the accumulated elastic energy G0 from tectonic loadings and the dissipation energy Gc due to frictional resistance. This theory predicts that rupture accelerates when G0 > Gc and decelerates when G0 < Gc, encapsulating a physics-based arrest criterion for dynamic ruptures: ∫ x x0+L (1 - Gc(x)/G0(x)) dx = 0, where L is predicted along-strike arrest location distance from the hypocenter (located at x0). The criterion applies to each of the rupture fronts on bilateral ruptures. For the AZF (Fig. 3) and DLSF (Fig. S3), the distribution of G0(x) along strike can be derived based on the inverted coupling map and the elapsed time since the last large earthquake and Gc(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 M6.8 earthquake, by assuming various hypocenter locations (Fig. 3B). We find that before the 1952 event, G0 on the asperity AZF-3 is higher than Gc (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 M6.8 earthquake in 1952 (Fig. 3B). After accounting for the stress perturbation of the 1952 earthquake (“Methods” section), we assess the G0, Gc 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 M6.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 M7.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 pressurization38 caused by the shear heating of fault zone pore fluids. The exponential friction law is also proposed in laboratory experiments39. 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 τr. Due to the lack of direct constraints, we need to indirectly infer τr value of the last big event. For the exponential friction, we assume τr equals the minimum dynamic strength τd (“Methods” section). Note that the absolute value of τd is irrelevant in elastic rupture dynamics on vertical planar faults9. For the power-law friction, we assume τr is larger than τd, corresponding to an empirical slip of the last big event (“Methods” section). As supported by laboratory experiments of large slips39, the assumption of τr = τ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 τr and thus result in 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 Gc scaling. 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 earthquakes40. 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 aftershocks42, 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 investigations43. 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/sb, where S and sb 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, G0 ≈ Gc. Because G0 increases faster than Gc 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 studies43,44. 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 studies26–28, 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.45 and combine them with the existing GNSS velocities46 for the kinematic inversion following the method of Wang and Shen46. 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 used26,

Vobs = VB + Vε + Vp + Vbackslip (1)

where Vobs are the observed GNSS velocities, VB and Vε the velocities due to steady state block rotation and interior strain26, respectively, Vp the postseismic viscoelastic relaxation effect of large historical earthquakes, and Vbackslip the deformation caused by fault coupling. VP 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 code47 with a layered infinite half-space model that consist 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 (VP 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) + Vis(x,t)H(t)] for slip deficit inversions (Eq. (2)), where H(t) is the Heaviside function, E(x) the elastic displacement, and 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 → ∞) using the layered viscoelastic model (Fig. S5).

VB and Vε are calculated by using an Euler vector and a strain rate tensor54 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 technique55. 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/a31. 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 VB, Vε, and VP from the observations, the residual deformation is assumed to be caused by the fault coupling and represented by Vbackslip. As in previous studies26–28,57, 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 method45. To make the solution stable, we add an a-priori smoothing constraint in the cost function,

F(sb) = ||Gsb - y||2 + α2||Hsb||2 (2)

where sb 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 control 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 = sb/s, where s is the fault slip rate. We estimate the uncertainty of the inverted slip deficit rate and coupling ratio using the statistic method53. 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” model54 and associated fault slip rate58,59. 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 estimate the b values homogeneously over space by employing the maximum-likelihood method within the ZMAP software61, 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 test62. For a good fit to the G-R law61, we only calculate the b values on the grids where the number of events (magnitude > Mc) exceeds Nmin. 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 tested the b values with different Nmin values (Fig. S10 and S11). Then, we further increase Nmin 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 Nmin values resulting in large uncertainties and high Nmin values leading to a highly smoothed variation (Figs. S10, S11). In this study, we therefore choose Nmin = 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 the fracture energy Gc, we compile the data from the literature63–66 and complement the data as inferred in natural faults38,67,68 and also observed in laboratory experiments69. As a first-order approximation derived from off-fault inelastic dissipation70,71 and thermal pressurization38, the scaling relation between the fracture energy and final slip can be expressed as a power-law, Gc ≈ BDn, where B and n are coefficients to be determined. Alternatively, laboratory experiments39 show that the friction weakening behaviors of natural rocks can be described by an exponential slip-weakening friction, τ = (τs - τd)e-D/dc + τd, where (τs - τd) and dc are the strength drop and the critical slip-weakening distance. For this friction, the scaling relation between fracture energy and final slip is Gc = (τs - τd)dc * (1 - (1 + D/dc)e-D/dc), where (τs - τd) and dc 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 τs - τd = 8 MPa and dc = 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 Gc constrained by dynamic models in the literature63–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/m3 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 zone48, the specific choice of the velocity model has a secondary effect on the numerical simulations. We use the spectral element software SPECFEM3D72 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: τ = τd + (τs - τd) / (1 + √(D/dc))p, where p is the power-law coefficient, τs, τd, and dc are the static strength, the dynamic strength, and the critical slip-weakening distance, respectively. In the simulations, we set p = 1 - n = 0.15, dc = 0.01 m, and τs - τd = 28 MPa, 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 τini = Δτaccu + τ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 τr = τd + (τs - τd) / (1 + √(5/dc))p.

For the exponential slip-weakening friction law, we assume τ = (τs - τd)e-D/dc + τd, where τs and τd, and dc 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 τs - τd = 8 MPa and dc = 0.2 m. Note that the absolute value of τd is not important in these two friction laws on vertical planer faults70, therefore here we artificially set τd = 1.2 MPa. In the exponential friction law, the initial shear stress is set as τini = Δτaccu + τd, where we assume that the latest large earthquakes dropped the shear stresses to the minimum dynamic strength τ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 Δτaccu, we first use the coupling ratio to calculate the shear stress changes per year τ and then calculate the total accumulated shear stress based on the elapsed time T since the latest large earthquakes, Δτaccu = τ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 width9,10. The 3D theory9,10 suggests a simple criterion for determining the arrest location of dynamic ruptures by evaluating the rupture potential: ϕ = ∫ x x0+L (1 - Gc(x)/G0(x)) dx = 0, where L is predicted along-strike arrest location distance from the hypocenter (located at x0), G0(x) and Gc(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 G0(x) is calculated based on the seismic coupling map and the elapsed time since the last M~7.5 large earthquakes. To obtain G0(x), we first estimate the total slip deficit Daccu based on the elapsed time T since the latest large earthquakes, and then calculate the accumulated shear stress Δτaccu on fault. The density of the potential elastic energy accumulated on the seismogenic fault is defined by 0.5Δτaccu Daccu. Following previous method10, we average 0.5Δτaccu Daccu across the seismogenic depth (i.e., W = 20 km in this paper) to obtain G0(x) along strike (Fig. 3A, C). The fracture energy Gc is estimated based on potential earthquake slip D and the scaling relation Gc(x) = BDn 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 use the accumulated slip deficit from the 1952 earthquake to the current time to calculate the potential fracture energy. We then calculate the density of the potential fracture energy by using BDnaccu and average it across the seismogenic depth to obtain Gc(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.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "General Comments:\nThe authors use newly-deployed GNSS stations that are located close to the fault traces of\nAnninghe and Zemuhe faults, and constructed a coupling model that is then validated with a\nvariety of external datasets. Then this model is used to evaluate the stress status along the\nfaults and are further used to assess the earthquake potential. Overall, the research is\ncomprehensive and results plus arguments are valid, though there are certain parts in the\nmanuscript that can be made better or clearer so that the audiences can better follow the\nstudy. Thus, I recommend a minor revision before it can be considered for publication.\nPlease see my detailed comments below.\nSome Detailed Comments:\nAbstract: what physics is used in the model? It is mentioned several times that the study\nused physics-based model, but it is not until the Methods section that is shown that the\nincluded physics is the viscoelastic effect. I suggest the audiences include some description\nof that in the early part of the manuscript. The audiences may feel confused.\nPlease use GNSS instead of GPS through-out the manuscript.\nI think the authors need to emphasize that the fault orientation makes it difficult for other\ngeodetic methods such as InSAR to detect the related fault motion. That’s also the reason\nwhy new GNSS occupations, especially considering that they are deployed right along the\nfaults, are quite valuable over this area.\nLine 65: 24 new-developed near-field GPS stations -> 24 newly-deployed near-field GNSS\nstations.\nFigure 1: Please plot the velocity vectors of the stations as well in 1a, and also mark the\nfaults in the inset figure too. Mark the location of historic events on the coupling map 1b, so\nas to show the spatial correlation with the inferred fault asperities. Plus, not all faults are\nproperly stated in the caption, e.g., XJHF. Make sure all of them are included in the caption.\nLine 68 and previous: The manuscript mentioned several times that not considering\nviscoelastic relaxation will result in overestimated fault coupling. Rather than just showing\ntwo models, I think the authors should include some discussion why viscoelastic effect may\nreduce coupling on fault. Do note NC is for broader audiences not just expert in geosciences.\nHow are the errorbars in figure 2a calculated? Also please explain what resulted in the\nnegative values in figure 2b, and how to digest these negative numbers.\nHow is high and low stress defined in figure 2b? The second from left panel seems to be not\nso low compared to the moderate panel on the right end.\nBased on the seismicity plotted on in figure 1b, I wonder whether there are enough of them\nto provide robust b value estimates. Could the authors add a figure in the supplements that\nshows how the b value is well estimated with low number of seismicity? e.g., N<100.\nLine 97: It is mentioned here that these low coupling areas might be permanent barriers for\nrupture propagation. Are there any analysis that support this argument. As I understand,\nlow coupling means the faults are less attached together, likely due to its location (at a low\nstress region) on a rough fault surface. In that sense, this part might be creeping or just\nmoves with rest part of the fault (stress shadow region). In that sense, it won’t be a barrier\non rupture propagation. I could be wrong on this but please elaborate a bit more in the\nmanuscript.\nThe simulation used in the model requires known initial stress status while the estimates is\nbased on the stress rate from the coupling model and the elapsed time since last large\nevent. I guess the underlying assumption here is that every past event has fully released the\naccumulated tectonic stress within its rupture range. This is a pretty strong assumption that\nI think the authors needs to make really clear in the manuscript and discuss the potential\nbias it may cause in the case the stress release is incomplete.\nI suggest that the authors use city+number to call the asperities or just use numbers to mark\nthem. Reading latitude in the manuscript and try to locate every asperity in the map is\npainful. You may use M-1 or M-2 for the two closes asperities to the city Mianing or X-1 and\nX-2 for the two ones closest to Xichang. Or alternatively, you may just use AZF-1, AZF-2 etc.\nAlso please mark these names in figure 1b.\nLine 187: when were the stations built?\nThe comparison shown in figure s5 is a little confusing. The difference between the two\nmodels is not just by adding XJHF, but also the block definition and data used seem different\ntoo. Please explain a bit in the methods section how these choices are made.\nDoes the block modeling used the new-deployed GNSS data? Maybe mark the with a\ndifferent color in figure S5.\nLine 261: There are multiple seismic velocity models over this region. Are these numbers\npicked based on any of the models?\nCode and data availability statement: It is recommended that the authors put their code on\na public accessible repository rather than making it available only upon request. I think\nNature publications requires the authors make materials, data, code and associated\nprotocols promptly available to readers without undue qualifications. This may not be\noptional.", "reviewer_2": "Dear Editor and authors,\nDiao et al invert fault coupling on the Anninghe-Zemuhe fault in southwestern China using\nGPS records from newly built and existing stations. They estimate rupture extents based on\nthe rupture potential calculated by the fracture energy Gc and potential elastic energy G0.\nG0 is calculated from the fault coupling and elapsed time since last major earthquakes. They\nalso simulate dynamic ruptures, calculate b values from recent seismicity to infer stress\nstates, to supplement their estimations on potential future earthquake magnitudes. I\nappreciate the simplicity of combining fault coupling ratio and elastic fracture theory to\nestimate rupture extents. The 24 new near-fault GPS stations and newly acquired records\nwould be valuable assets to the earthquake community. The paleoseismic data suggest a\nmixture of single- and multi-segment (in terms of coupling ratio) ruptures on the Anninghe-\nZemuhe fault in the past 500 years. However, the static nature of the rupture potential\nmethod may not be able to evaluate multi-segment ruptures. Dynamic rupture model is\nable to do so but the planar fault assumption may not properly represent the fault\ngeometric complexities in the fault system. It is also surprising to see that fault geometric\ncomplexities, which are barriers to ruptures given geologic and numerical evidence, are not\ntouched in the manuscript. The new data and inverted fault coupling are definitely valuable,\nbut I am afraid the rupture potential method may provide incremental scientific merit in\nevaluating rupture extents and their likelihoods in complex fault networks compared to\nestablished literature.\nMy major and detailed minor comments are provided below.\nMajor comments:\n1. Fault geometric complexities such as fault bends and step-overs are associated with\nrupture terminations and initiations (e.g., Wesnousky, 2006, Biasi and Wesnousky, 2017).\nProgress has been made in the last few decades by combing physics-based dynamic rupture\nand earthquake cycle models to address conditions and likelihood for multi-segment\nruptures, which I think is the focus of this study too. However, the introduction seems to\ncover few of them. For example, dynamic rupture model is a powerful tool to explore\nconditions for multi-segment ruptures, which could greatly increase the sizes (e.g., Harris et\nal., 1991, Harris and Day, 1993, Ulrich et al., 2019). Dynamic ruptures constrained by\ngeodetically inverted coupling provides insights on rupture segmentation and potential\nscenarios for seismic hazard assessment (e.g., Ramos et al., 2021 and the cited Yang et al.,\n2019). Going beyond single dynamic ruptures, earthquake cycle models with fault geometric\ncomplexity demonstrate that ruptures conditionally break geometric barriers, leading to\nvariable earthquake sizes (e.g., Duan and Oglesby, 2005). With parameterized dynamic\nrupture effect, physics-based earthquake simulator such as RSQsim could produce\nearthquake sequences of various sizes and reproduce seismic statistics (e.g., Dieterich and\nRichards-Dinger, 2010). Applied to complex fault systems in Southern California, the\nsimulator improves the Probabilistic Seismic Hazard Analysis (PHSA) rupture models,\nallowing physics-based estimations of multi-segment ruptures (Shaw et al., 2018).\nIt looks like the strike of ANH-ZMH fault changes much and DLS fault is quite segmented. I\nwonder how the fault geometric features would contribute to the understanding of rupture\nsegmentation in this fault system, and complement the findings from this work.\n2. The rupture potential method is static. I am not sure if the static method could address\nthe occurrence of multi-segment M 7.5 ruptures shown in the paleoseismic datasets. In\nFigure 2a, ruptures that connect two high coupling segments (1480 M7.5, 1536 M7.5, 1850 M7.5) are mixed with ruptures limited within single high coupling segment (1732 M6.8, 1952 M6.8, and potentially 1489 M6.8). There are also interesting dynamics that the 1489 M 6.8 was followed shortly be the 1536 M7.5 rupture, and that the 1732 M 6.8 rupture was\nfollowed by the 1850 M 7.5 rupture. Dynamic rupture model is an adequate tool to study\nmulti-segment ruptures. However, if I understand the dynamic rupture models correctly, a\nplanar vertical fault is used, which seems to misrepresent the Anninghe-Zemuhe fault\nsystem with varying strike.\n3. More details on the model setup, choices of parameters, equations of the methods would\nbe helpful. Details are scattered across the manuscript, figure captions, and supplementary\nmaterial. For example, how is the fault trace determined? How wide is the fault? Is the block\nmodel used in inversions 3-D and share similar velocity structures with the viscoelastic\nmodel?\nMore detailed comments:\nLine 23: I think the numerical model may supplement the theoretical model but not validate\nit. You need observations to validate a model.\nLine 37-38. It seems the supercycle concept is mostly used for subduction earthquakes but\nnot much applied to complex networks of continental faults such as those in Southern\nCalifornia.\nLine 41. And uncertainties in stress states?\nLine 71. What is the geologic estimation of slip rate on DLSF? And how does the geodetic\nslip rate compare to observations?\nLine 77-86. Not sure if earthquake sizes can be inferred from fault maturities. It may be\nnecessary to compare fault systems with similar geologic accumulated slip, age, and tectonic\nloading condition, and earthquake histories to make such calls.\nLine 127. Are you simulating the whole Anninghe-Zemuhe fault in dynamic rupture models?\nIf so, how do you justify the planar fault assumption given the change in strike.\nLine 143-144. It would be nice to show rupture extents from the models and compare them\nto the coupling ratio distribution. Are there any multi-segment ruptures comparable to\nthose paleoseismic M 7.5 events? Are there even larger events allowed by the models?\nLine 188: the GPS records are critical to the inversion. It would be helpful to talk about, even\nbriefly, uncertainties and qualities of GPS records. For newly deployed stations, it would be\nhelpful to have a table or data file to show the station locations, EW and NS rates and\nuncertainties for the new data. Are there any vertical components?\nFor Figure S5, is this block model 2-D or 3-D? Block models are known to suffer trade-offs\nbetween fault widths and inverted slip rate. Is there such an issue in this study?\nFigure S5, would be nice to have fault names on the maps.\nLine 200. While Model 1 that ignores the XJHF generates a slip rate on AZF close to geologic\nobservations, is there any estimation of geologic slip rate on the XJHF? If so, how does it\ncompare to AZF? Given the close proximity of the two faults, can XJHF be safely ignored?\nLine 201. It would be helpful to briefly introduce the geologic slip rate estimations and\nuncertainties for the faults involved.\nLine 203. What is the model setup of the viscoelastic model in terms of model domain, fault\ngeometry and historical earthquake slip implementation? What about the availability of\ncode PSGRN/PSCMP and the viscoelastic model?\nLine 204. Do you mean Figure S6?\nLine 212. Is s backslip or backslip rate?\nLine 214. Do you have the trade-off plot between \\alpha and data misfit?\nLine 215. Would be helpful to lay out equations used to calculate fault coupling ratio from\ninverted slip rates. Also, would be helpful to show briefly how the uncertainties are\nevaluated.\nLine 216-224. Seems this paragraph that discusses viscoelastic effect could be merged with\nLine 202-205?\nLine 257. It seems a planar fault geometry is used. Is the whole Anninghe-Zemuhe fault\nsystem simulated in dynamic ruptures? The fault strike variations are not small to be\nignored. What is the total faut length?\nLine 262. What is the CFL coefficient and time step size? Would be helpful to have a table of\nkey parameters.\nLine 265. Same to the comment to Line 23.\nData Availability. Is there any public data repository to store and share the data?\nCode availability. Nature (and other journals in general) recommends DOI-minting\nrepository to store and share code. What about the code availability of PSGRN/PSCMP? Is\nComputational Infrastructure a data repository?\nBiasi, G.P. & Wesnousky, S.G., 2017. Bends and Ends of Surface Ruptures, Bulletin of the\nSeismological Society of America, 107, 2543-2560.\nDieterich, J.H. & Richards-Dinger, K.B., 2010. Earthquake Recurrence in Simulated Fault\nSystems, Pure and Applied Geophysics, 167, 1087-1104.\nDuan, B. & Oglesby, D.D., 2005. Multicycle dynamics of nonplanar strike‐slip faults, Journal\nof Geophysical Research: Solid Earth (1978–2012), 110.\nHarris, R.A., Archuleta, R.J. & Day, S.M., 1991. Fault steps and the dynamic rupture process:\n2-D numerical simulations of a spontaneously propagating shear fracture, Geophysical\nResearch Letters, 18, 893-896.\nHarris, R.A. & Day, S.M., 1993. Dynamics of fault interaction: parallel strike-slip faults,\nJournal of Geophysical Research, 98, 4461.\nRamos, M.D., Huang, Y., Ulrich, T., Li, D., Gabriel, A.-A. & Thomas, A.M., 2021. Assessing\nMargin-Wide Rupture Behaviors Along the Cascadia Megathrust With 3-D Dynamic Rupture\nSimulations, Journal of Geophysical Research: Solid Earth, 126, e2021JB022005.\nShaw, B.E., Milner, K.R., Field, E.H., Richards-Dinger, K., Gilchrist, J.J., Dieterich, J.H. &\nJordan, T.H., 2018. A physics-based earthquake simulator replicates seismic hazard statistics\nacross California, 4, eaau0688.\nUlrich, T., Gabriel, A.A., Ampuero, J.P. & Xu, W., 2019. Dynamic viability of the 2016 Mw 7.8\nKaikoura earthquake cascade on weak crustal faults, Nat Commun, 10, 1213.\nWesnousky, S.G., 2006. Predicting the endpoints of earthquake ruptures, Nature, 444, 358-\n360.", "reviewer_2_round_2": "I appreciate the authors for their diligent efforts in addressing my comments.\nThe manuscript is much improved. I endorse its publication now. Thanks!\nOne last comment: I still feel there is room for improvement related to my previous\ncomment \"In Figure 2a, ruptures that connect two high coupling segments (1480 M7.5, 1536 M7.5, 1850 M7.5)\nare mixed with ruptures limited within single high coupling segment (1732 M6.8, 1952 M6.8, and\npotentially 1489 M6.8). There are also interesting dynamics that the 1489 M 6.8 was followed shortly\nbe the 1536 M7.5 rupture, and that the 1732 M 6.8 rupture was followed by the 1850 M 7.5\nrupture.\nDynamic rupture model is an adequate tool to study multi-segment ruptures.\"\nThese potential multi-segment ruptures are relevant to seismic hazard evaluation. I look\nforward to seeing progress from your method to address the potentials of multi-segment\nruptures, and the complex dynamics in terms of both their rupture extents and recurrences.\nFault geometric complexity affects rupture extents through heterogeneous stress conditions\nthat are inherited from past earthquakes and developed over cycles. In Figure 4C, it appears\nthat ruptures are well confined within bend limits. So L184 may not be true. Anyhow, they\ndon't contradict to your method, because the heterogeneous stress states are well\ncharacterized through your inversion with GPS constraints. It is nice to see a planar fault\nperforming so well with your method."}, "author_response": "Response to reviewers' comments: Black Italic = Comment, Blue = Response,\nRed = Revision\n===================================================================\nReviewer #1 (Remarks to the Author):\n1. The authors use newly-deployed GNSS stations that are located close to the fault traces of Anninghe\nand Zemuhe faults, and constructed a coupling model that is then validated with a variety of external\ndatasets. Then this model is used to evaluate the stress status along the faults and are further used to\nassess the earthquake potential. Overall, the research is comprehensive and results plus arguments are\nvalid, though there are certain parts in the manuscript that can be made better or clearer so that the\naudiences can better follow the study. Thus, I recommend a minor revision before it can be considered\nfor publication. Please see my detailed comments below.\nResponse: Thank you very much for the valuable comments and suggestions, which helped us improve\nthis manuscript. Please find in the following our point-to-point responses to the specific comments.\n2. Abstract: what physics is used in the model? It is mentioned several times that the study used physics-\nbased model, but it is not until the Methods section that is shown that the included physics is the\nviscoelastic effect. I suggest the audiences include some description of that in the early part of the\nmanuscript. The audiences may feel confused.\nResponse: The term ‘physics’ refers to the 3D rupture dynamics of faults. To clarify it, we revised in the\nabstract:\nHere, we show how the maximum possible earthquake magnitude of a fault system can be assessed by\ncombining high-resolution fault coupling maps with a physics-based theoretical model derived from\nthree-dimensional dynamic fracture mechanics, as confirmed by dynamic rupture simulations.\nand in Introduction:\nHere, we combine kinematic inversions and a physics-based theoretical model based on three-dimension\n(3D) dynamic fracture mechanics to show that…\nand in the main text:\nThe accumulated slip deficit and shear stress, combined with 3D dynamic rupture theory, allow us\nto assess the maximum possible earthquake magnitude on the AZF and the DLSF. In dynamic rupture\ntheory9, rupture propagation and arrest on a fault are controlled by the accumulated elastic\nenergy 𝐺 from tectonic loadings and the dissipation energy 𝐺 due to frictional resistance. This theory\n0 𝑐\npredicts that rupture accelerates when 𝐺 > 𝐺 and decelerates when 𝐺 < 𝐺 , encapsulating a physics-\n0 𝑐 0 𝑐\nbased arrest criterion for dynamic ruptures: ∫\n𝑥0+𝐿\n(1−𝐺 (𝑥)/𝐺 (𝑥))𝑑𝑥 = 0 , where 𝐿 is predicted\n𝑥0 𝑐 0\nalong-strike arrest location distance from the hypocentre (located at 𝑥 ). The criterion applies to each\n0\nof the rupture fronts on bilateral ruptures. For the AZF (Fig. 3) and DLSF (Fig. S3), the distribution of\n𝐺 (𝑥)along strike can be derived based on the inverted coupling map and the elapsed time since the last\n0\nlarge earthquake and 𝐺 (𝑥) can be determined using a power-law slip-weakening friction model to fit\n𝑐\nglobal earthquake observations (Methods).\n3. Please use GNSS instead of GPS through-out the manuscript.\nResponse: We have used ‘GNSS’ throughout the manuscript.\n4. I think the authors need to emphasize that the fault orientation makes it difficult for other geodetic\nmethods such as InSAR to detect the related fault motion. That’s also the reason why new GNSS\noccupations, especially considering that they are deployed right along the faults, are quite valuable over\nthis area.\nResponse: Thanks for this suggestion. We added the following sentence in the revised main text:\nBecause the AZF has a near N-S trending that limits the applicability of the InSAR\n(Interferometric Synthetic Aperture Radar) technique, the near-field GNSS observations in this study can\nprovide valuable constraints on fault coupling distribution.\n5. Line 65: 24 new-developed near-field GPS stations -> 24 newly-deployed near-field GNSS stations.\nResponse: Revised.\n6. Figure 1: Please plot the velocity vectors of the stations as well in 1a, and also mark the\nfaults in the inset figure too. Mark the location of historic events on the coupling map 1b, so as to show the spatial\ncorrelation with the inferred fault asperities. Plus, not all faults are properly stated in the caption,\ne.g., XJHF. Make sure all of them are included in the caption.\nResponse: Thanks for these comments. We have revised Fig. 1 and the figure caption accordingly:\nFigure 1. Tectonic setting and fault coupling ratio on the Anninghe-Zemuhe fault system. A.\nTectonic setting of the study area. Black solid lines indicate faults. The red and gray triangles represent\nthe 24 new-developed near-field GNSS stations and existing stations used for constraining the kinematic\ncoupling model. Blue stars indicate historical events (M>6.5) in the past 500 years. Blue vectors show\nthe GNSS velocity in this region. ANHF: Anninghe fault; ZMHF: Zemuhe fault; DLSF: Daliangshan\nfault; LJ-XJHF: Lijiang-Xiaojinhe fault; LMSF: Longmen Shan fault. B. Coupling ratio on the AZF and\nDLSF. The fault changes from a completely locked state to a fully steady slip state as the coupling ratio\nvaries from 1 to 0. Gray dots show seismicity distribution that was projected onto the faults. AZF-1 to\nAZF-6 show the identified asperities on the AZF.\n7. Line 68 and previous: The manuscript mentioned several times that not considering\nviscoelastic relaxation will result in overestimated fault coupling. Rather than just showing\ntwo models, I think the authors should include some discussion why viscoelastic effect may\nreduce coupling on fault. Do note NC is for broader audiences not just expert in geosciences.\nResponse: Thanks for the comment. We explained the effect of viscoelastic relaxation on fault coupling\ninversions in the revised main text：\nIn the interseismic stage, viscoelastic relaxation at depth can affect crustal deformations in middle- and\nfar-filed areas29. Ignoring this effect in the kinematic inversions results in larger slip deficits on deep\nfault segments to compensate, which thus leads to an overestimation of their coupling ratio30. Therefore,\nwe adopt viscoelastic deformation models to infer the coupling ratio of the AZF and DLSF. Our results\nfurther confirm that the elastic model overestimates the fault coupling on deep patches (Fig. S1).\n8. How are the errorbars in figure 2a calculated?\nResponse: We added explanation in Methods:\nWe estimate the uncertainty of the inverted slip deficit rate and coupling ratio using the statistic method53.\nThe details are as follows: we first generate a set of Gaussian noise signals using the observation\nuncertainty as the standard deviation. Then, we add them to the observation data to get several sets of\ndisturbed synthetic data. For each of the disturbed data sets, we repeated the inversion process to get the\ncoupling ratio on each fault patch and the mean coupling ratio above 15 km. These steps were run 500\ntimes in order to estimate the uncertainty of the parameters statistically.\nAlso please explain what resulted in the negative values in figure 2b, and how to digest these negative\nnumbers.\nResponse: The negative values in Figure 2B, primarily near the bottom of the asperities, are caused by\nlarge slip gradient near the edge of locked asperities. These could be mitigated by using denser patches\nand stronger spatial smoothing, nevertheless it requires additional constraints and assumptions for the\ndeep patches. However, we believe that the effects of these negative values are relatively minor for\nseismic hazard assessment, as the dense near-fault GNSS observations can provide very good first-order\nconstraints on seismogenic fault zones (see the resolution tests shown in Fig. S2).\n9. How is high and low stress defined in figure 2b? The second from left panel seems to be not\nso low compared to the moderate panel on the right end.\nResponse: These stresses are defined qualitatively. To avoid confusion, we have revised figure 2b by\nlabeling the specific values of the stresses.\n10. Based on the seismicity plotted on in figure 1b, I wonder whether there are enough of them\nto provide robust b value estimates. Could the authors add a figure in the supplements that\nshows how the b value is well estimated with low number of seismicity? e.g., N<100.\nResponse: The estimates of the b value in figure 1b are robust. To verify this, we conducted two\nadditional tests, shown in Figures S10 and S11, and added new explanations in the main text:\nTo verify the robustness of the inferred b values, we conducted additional tests. First, we increase 𝑁\n𝑚𝑖𝑛\nfrom 30 to 60 and observe the resulting variation in the calculated b values (Fig. S10). Then, we further\nincrease 𝑁\n𝑚𝑖𝑛\nto 80 and 100, but with an enlarged sampling radius r of 40 km (Fig. S11). We find that\nthe primary features of the inferred b values are qualitatively consistent, despite low 𝑁 values\n𝑚𝑖𝑛\nresulting in large uncertainties and high 𝑁 values leading to a highly smoothed variation (Fig. S10\n𝑚𝑖𝑛\n& S11). In this study, we choose 𝑁 =50 with a sampling radius r of 25 km, following ref37, to ensure\n𝑚𝑖𝑛\nsufficient spatial sampling while preserving spatial variations of the b values in each asperity. In addition,\nthe frequency–magnitude distributions within selected patches also show that the estimates of the b\nvalue are robust (Fig. S10).\nFigure S10. A-D show the variation in the calculated b values along the Anninghe-Zemuhe fault for\ndifferent 𝑁 that represents the defined minimum number of events with magnitude above Mc. E\n𝑚𝑖𝑛\nshows the frequency–magnitude distribution for annotated patches indicated by 1-4 in Fig. C.\nFigure S11. The b values along the Anninghe-Zemuhe fault calculated with variable r and 𝑁 (in\n𝑚𝑖𝑛\nthe legend).\n11. Line 97: It is mentioned here that these low coupling areas might be permanent barriers for\nrupture propagation. Are there any analysis that support this argument. As I understand,\nlow coupling means the faults are less attached together, likely due to its location (at a low\nstress region) on a rough fault surface. In that sense, this part might be creeping or just\nmoves with rest part of the fault (stress shadow region). In that sense, it won’t be a barrier\non rupture propagation. I could be wrong on this but please elaborate a bit more in the\nmanuscript.\nResponse: Thanks for point this out. The term \"permanent barrier\" usually indicates a region that never\nslips seismically (i.e., only slips aseismically). However, we used double quotation marks (\"permanent\")\nbecause we believe these barriers are long-standing rather than truly permanent. To avoid unexpected\nconfusion, we have replaced “permanent” with “long-standing” in that sentence.\n12. The simulation used in the model requires known initial stress status while the estimates is\nbased on the stress rate from the coupling model and the elapsed time since last large\nevent. I guess the underlying assumption here is that every past event has fully released the\naccumulated tectonic stress within its rupture range. This is a pretty strong assumption that\nI think the authors needs to make really clear in the manuscript and discuss the potential\nbias it may cause in the case the stress release is incomplete.\nResponse: To clarify this point, we added in the main text:\nIn the numerical model, we implement the power and exponential slip-weakening friction laws,\nrespectively, whose coefficients are determined based on the data of global earthquakes within the\nrelevant range of slip values (Methods). These two friction laws represent two end-member regimes of\nthermal pressurization38 caused by shear heating of fault zone pore fluids. The exponential friction law\nis also proposed in laboratory experiments39. The current shear stress on fault has been building up since\nthe latest large earthquake, during which the frictional strength dropped to a residual value 𝜏 .\n𝑟\nDue to lack of direct constraints, we need to indirectly infer 𝜏 value of the last big event. For the exponential\n𝑟\nfriction, we assume 𝜏 equals the minimum dynamic strength 𝜏 (Methods). Note that the absolute\n𝑟 𝑑\nvalue of 𝜏 is irrelevant in elastic rupture dynamics on vertical planar faults9. For the power-law friction,\n𝑑\nwe assume 𝜏 is larger than 𝜏 , corresponding to an empirical slip of the last big event (Methods). As\n𝑟 𝑑\nsupported by laboratory experiments of large slips39, the assumption of 𝜏 = 𝜏 for the exponential\n𝑟 𝑑\nfriction relies on the estimated slip of the last big event being sufficiently large. The assumption for the\npower friction relies on the estimated slip of the historical big event: overestimated slips lead to\nunderestimated 𝜏 and thus result in underestimated potential earthquake sizes. Despite these different\n𝑟\nassumptions of shear stresses and friction laws, we find that theoretical and numerical predictions (Fig.\n3B and 3D) are quantitatively consistent in different rupture scenarios within the uncertainties in fitting\nthe G scaling data.\n𝑐\n13. I suggest that the authors use city+number to call the asperities or just use numbers to mark\nthem. Reading latitude in the manuscript and try to locate every asperity in the map is\npainful. You may use M-1 or M-2 for the two closes asperities to the city Mianing or X-1 and\nX-2 for the two ones closest to Xichang. Or alternatively, you may just use AZF-1, AZF-2 etc.\nAlso please mark these names in figure 1b.\nResponse: Revised following the comment. We use AZF-1, AZF-2, etc., thorough this paper.\n14. Line 187: when were the stations built?\nResponse: We add the information in the main text：\nTo better constrain the coupling ratio distribution along the AZF, we built 24 new near-fault GNSS\nstations in 2018 and surveyed them three times in 2019, 2020 and 2022 (Fig. 1A).\n15. The comparison shown in figure s5 is a little confusing. The difference between the two\nmodels is not just by adding XJHF, but also the block definition and data used seem different\ntoo. Please explain a bit in the methods section how these choices are made.\nResponse: Thanks for the comment. We explained in Methods:\nWe divide the study area into three blocks (Model 1 in Fig. S7) based on the strategy inferred from a\nEuler pole clustering technique55. We also test the other block division strategy (Model 2 in Fig. S7) to\ninvestigate whether the LJ-XJHF shall be defined as the block boundary in the model. In Model 1 and\nModel 2, we use different coverages of GNSS stations to define the motion of blocks west of the AZF,\ni.e. Chuandian block and Dianzhong block in Fig. S7. Modeling results indicate that Model 2 leads to a\nslip rate of 9.2-5.7mm/a along the AZF (Model 2 in Fig. S7) that is higher than the geological result of\n6.5±1.0 mm/a22. On the contrary, Model 1 that neglects the XJHF, generates a slip rate of AZF (6.1-4.9\nmm/a) that agrees better with the geologic investigations.\nand in the caption of Fig. S7 (previous Fig. S5):\nBlock division strategy, inferred fault slip rates and model residual. Models 1 and 2 show block division\nstrategies that exclude and include the Xiaojinhe fault, respectively. Left: block division strategy used\nfor slip deficit inversion. Vectors within the regions outlined by black dashed lines were used to estimate\nthe parameter of block motions. Middle: model residuals and resulting slip rate distribution. Right:\nhistograms of residuals in EW and NS components. CD: Chuandian block; DLS: Daliangshan block;\nMB: Mabian block; DZ: Dianzhong block.\n16. Does the block modeling used the new-deployed GNSS data? Maybe mark the with a\ndifferent color in figure S5.\nResponse: The new-deployed GNSS data are close to the fault and thus are not used in the block model,\nbecause near-fault stations are mainly controlled by fault locking/creep. We clarified this point in the\ncaption of Fig. S7 (previous S5).\nVectors within the regions outlined by black dashed lines were used to estimate the parameter of block\nmotions.\n17. Line 261: There are multiple seismic velocity models over this region. Are these numbers\npicked based on any of the models?\nResponse: We explained in the method:\nThe P wave speed, S wave speed and density of the medium of the elastic model are uniformly assumed\nas 5770 m/s, 3330 m/s and 2705 g/m3 based on parameters of upper layers shown by Liu et al48. As the\nvariation of shear modulus is less than 10% within the seismogenic zone48, the specific choice of the\nvelocity models has a secondary effect on the numerical simulations.\n18. Code and data availability statement: It is recommended that the authors put their code on\na public accessible repository rather than making it available only upon request. I think\nNature publications requires the authors make materials, data, code and associated\nprotocols promptly available to readers without undue qualifications. This may not be\noptional.\nResponse: Thanks for the comments. We addressed this issue by adding the following information in\nthe ‘Data availability’\nThe velocity of newly-developed GNSS stations is listed in Table S1, and the other GNSS velocities can\nbe found in the supplementary files of ref 46. The refined seismic catalog used for b value calculation is\nshown in Supplementary Data 1.\nand ‘Code availability’:\nThe PSGRN/PSCMP code for postseismic deformation modeling and viscoelastic Green’s function\ncalculation, as well as the SDM code for the slip deficit inversions, are available at the public sever of\nGerman Research Centre for Geosciences (GFZ): ftp://ftp.gfz-potsdam.de/pub/home/turk/wang/. The\nopen-source software SPECFEM3D used in the 3D fully dynamic rupture simulations are available from\nthe Computational Infrastructure: https://github.com/SPECFEM/specfem3d. The open-source software\nZMAP used for b value calculations is available at: https://github.com/zmap/zmap.\nResponse to reviewers' comments: Black Italic = Comment, Blue = Response,\nRed = Revision\n===================================================================\nReviewer #2 (Remarks to the Author):\n0.Diao et al invert fault coupling on the Anninghe-Zemuhe fault in southwestern China using GPS\nrecords from newly built and existing stations. They estimate rupture extents based on the rupture\npotential calculated by the fracture energy Gc and potential elastic energy G0. G0 is calculated from\nthe fault coupling and elapsed time since last major earthquakes. They also simulate dynamic ruptures,\ncalculate b values from recent seismicity to infer stress states, to supplement their estimations on\npotential future earthquake magnitudes. I appreciate the simplicity of combining fault coupling ratio\nand elastic fracture theory to estimate rupture extents. The 24 new near-fault GPS stations and newly\nacquired records would be valuable assets to the earthquake community. The paleoseismic data suggest\na mixture of single- and multi-segment (in terms of coupling ratio) ruptures on the Anninghe-Zemuhe\nfault in the past 500 years. However, the static nature of the rupture potential method may not be able\nto evaluate multi-segment ruptures. Dynamic rupture model is able to do so but the planar fault\nassumption may not properly represent the fault geometric complexities in the fault system. It is also\nsurprising to see that fault geometric complexities, which are barriers to ruptures given geologic and\nnumerical evidence, are not touched in the manuscript. The new data and inverted fault coupling are\ndefinitely valuable, but I am afraid the rupture potential method may provide incremental scientific merit\nin evaluating rupture extents and their likelihoods in complex fault networks compared to established\nliterature.\nResponse: Thanks for the inspiring and constructive comments, which helped us improve this\nmanuscript. Regarding the primary concerns about the rupture potential method and the complexities of\nfault geometry, we summarize our responses to these two points as follows: 1) The rupture potential\nmethod is dynamic rather than static and is originally derived from the 3D theory of dynamic rupture\nmechanics; 2) Our additional 3D dynamic rupture simulations with curved faults, as shown in newly\nadded Figure 4, demonstrate that the geometrical effects are insignificant in predicting the maximum\npossible earthquake magnitudes on the AZF. The additional tests show that the rupture potential method\ncan serve as theoretical basis for evaluating rupture extents on a fault system with dense near-fault\ngeodetic data. In the revised version, we believe that all the comments and questions have been\nadequately addressed. Please find below point-to-point responses to each specific comments.\n1. Fault geometric complexities such as fault bends and step-overs are associated with rupture\nterminations and initiations (e.g., Wesnousky, 2006, Biasi and Wesnousky, 2017). Progress has been\nmade in the last few decades by combining physics-based dynamic rupture and earthquake cycle models\nto address conditions and likelihood for multi-segment ruptures, which I think is the focus of this study\ntoo. However, the introduction seems to cover few of them. For example, dynamic rupture model is a\npowerful tool to explore conditions for multi-segment ruptures, which could greatly increase the sizes\n(e.g., Harris et al., 1991, Harris and Day, 1993, Ulrich et al., 2019). Dynamic ruptures constrained by\ngeodetically inverted coupling provides insights on rupture segmentation and potential scenarios for\nseismic hazard assessment (e.g., Ramos et al., 2021 and the cited Yang et al., 2019). Going beyond single\ndynamic ruptures, earthquake cycle models with fault geometric complexity demonstrate that ruptures\nconditionally break geometric barriers, leading to variable earthquake sizes (e.g., Duan and Oglesby,\n2005). With parameterized dynamic rupture effect, physics-based earthquake simulator such as RSQsim\ncould produce earthquake sequences of various sizes and reproduce seismic statistics (e.g., Dieterich\nand Richards-Dinger, 2010). Applied to complex fault systems in Southern California, the simulator\nimproves the Probabilistic Seismic Hazard Analysis (PHSA) rupture models, allowing physics-based\nestimations of multi-segment ruptures (Shaw et al., 2018).\nResponse: Thank you for pointing this out. We have acknowledged these developments in Introduction:\nAn accessible way to assess future seismic hazard is to estimate the maximum possible earthquake\nmagnitudes on faults as a function of elapsed time, capturing for example the possibility that a fault\nmight not be ready to break over its whole length. Kinematically-based11-13 or dynamically-based8,14-20\nmodels on various fault systems, including estimates of their slip budget, have been used to infer their\npotential seismic moments. Specifically, fault geometric complexities21,22, such as fault bends and step-\novers, have been integrated into single rupture simulations15-18 and earthquake cycle models19,20.\nFurthermore, complex fault systems have also been incorporated with a physics-based statistical\nsimulator23 to reproduce seismic hazard statistics in Southern California. Both kinematically-based and\ndynamically-based models provide important insights on rupture segmentation and seismic hazard\nassessment. However, kinematically-based methods, although efficient, may overestimate the maximum\nmagnitudes of potential earthquakes due to lack of integration with rupture dynamics, while\ndynamically-based numerical simulations are too computationally expensive for efficient hazard\nassessment and thorough sensitivity analysis.\nIt looks like the strike of ANH-ZMH fault changes much and DLS fault is quite segmented. I wonder how\nthe fault geometric features would contribute to the understanding of rupture segmentation in this fault\nsystem, and complement the findings from this work.\nResponse: Thanks for this insightful comment. We have performed additional 3D dynamic rupture\nsimulations with curved faults, as shown in newly added Figure 4, and found that the geometrical effects\ndo not affect our conclusion. We added a new paragraph in the main text:\n...\nWe then conduct additional 3D dynamic rupture models with a curved fault as used in the coupling\ninversion model (red curve in Fig. 4A), characterized by a significant bending within segment AZF-4.\nFor comparison, we also consider a model with smoothed fault trace (blue curve in Fig. 4A). We find\nthat fault bending is not the primary factor controlling earthquake sizes on the AZF, and that the\nnumerical simulations are quantitatively consistent with the theoretical prediction regardless of different\nfault geometries on the AZF (Fig. 4B & 4C). The explanation may be that the dominant wavelength of\nfault curvature on the AZF is larger than the length scales of rupture depths or seismic asperities.\nConsequently, the first-order rupture dynamics on the AZF can be well predicted by 3D dynamic fracture\nmechanics based on a planar fault. While in the future, it is necessary to consider the effects of fault\ncomplexities, such as step-overs on the DLSF, despite its observations being less constrained compared\nto those of the AZF.\nFigure 4. Effects of fault geometries on maximum possible magnitude of earthquakes on the AZF.\nA. The surface traces of curved and planar faults (legend) used in dynamic rupture simulations. B. The\nmaximum possible earthquake magnitude in 2023 as a function of hypocenter location predicted by\nsimulations with curved and planar faults (legend in A). The simulation results are compared with the\ntheoretical prediction. C. Accumulated shear stress on curved fault (top) and superimposition of slip for\nfour earthquake scenarios on curved (upper middle), smoothed (lower middle), and planar (bottom) fault,\neach assuming a different nucleation location (marked by purple stars).\n2. The rupture potential method is static. I am not sure if the static method could address the occurrence\nof multi-segment M 7.5 ruptures shown in the paleoseismic datasets.\nResponse: The rupture potential method is dynamic rather than static. To clarify this point, we added in\nthe main text:\nThe rupture potential method is dynamic rather than static. To clarify this point, we added in the main text:\nThe accumulated slip deficit and shear stress, combined with 3D dynamic rupture theory,\nallows us to assess the maximum possible earthquake magnitude on the AZF and DLSF. In dynamic\nrupture theory9, rupture propagation and arrest on a fault are controlled by the accumulated\nelastic energy 𝐺 from tectonic loadings and the dissipation energy 𝐺 due to frictional\n0 𝑐\nresistance. This theory predicts that rupture accelerates when 𝐺 > 𝐺 and decelerates\n0 𝑐\nwhen 𝐺 < 𝐺 , encapsulating a physics-based arrest criterion for dynamic ruptures: ∫\n𝑥0+𝐿\n(1−\n0 𝑐 𝑥0\n𝐺 (𝑥)/𝐺 (𝑥))𝑑𝑥 = 0, where 𝐿 is predicted along-strike arrest location distance from the\n𝑐 0\nhypocentre (located at 𝑥 ). The criterion applies to each of the rupture fronts on bilateral\n0\nruptures. For the AZF (Fig. 3) and DLSF (Fig. S3), the distribution of 𝐺 (𝑥) along strike can be\n0\nderived based on the inverted coupling map and the elapsed time since the last large earthquake\nand 𝐺 (𝑥) can be determined using a power-law slip-weakening friction model to fit global\n𝑐\nearthquake observations (Methods).\nIn Figure 2a, ruptures that connect two high coupling segments (1480 M7.5, 1536 M7.5, 1850 M7.5)\nare mixed with ruptures limited within single high coupling segment (1732 M6.8, 1952 M6.8, and\npotentially 1489 M6.8). There are also interesting dynamics that the 1489 M 6.8 was followed shortly\nbe the 1536 M7.5 rupture, and that the 1732 M 6.8 rupture was followed by the 1850 M 7.5 rupture.\nDynamic rupture model is an adequate tool to study multi-segment ruptures. However, if I understand\nthe dynamic rupture models correctly, a planar vertical fault is used, which seems to misrepresent the\nAnninghe-Zemuhe fault system with varying strike.\nResponse: Please refer to our response to comment 1.\n3. More details on the model setup, choices of parameters, equations of the methods would be helpful.\nDetails are scattered across the manuscript, figure captions, and supplementary material.\nResponse: We have added more details of both kinematic and dynamic models. For kinematic models,\nplease refer to the responses to comments 11, 12, and 16. For dynamic models, please refer to the\nresponses to comments 23.\nFor example, how is the fault trace determined?\nResponse: We clarified this in the main text:\nWe set the fault trace based on the geologic data56 and assume the AZF to be vertical following previous\nstudies26-28, 57.\nHow wide is the fault?\nResponse: The width of the seismogenic fault is 20 km.\nIs the block model used in inversions 3-D and share similar velocity structures with the viscoelastic\nmodel?\nResponse: Yes. They are the same.\n4. Line 23: I think the numerical model may supplement the theoretical model but not validate it. You\nneed observations to validate a model.\nResponse: We rephrased it by \"confirm\".\n5. Line 37-38. It seems the supercycle concept is mostly used for subduction earthquakes but not much\napplied to complex networks of continental faults such as those in Southern California.\nResponse: Thanks for point this out. We revised it as “subduction zone earthquakes”.\n6. Line 41. And uncertainties in stress states?\nResponse: Thanks, we added in the main text:\nPredicting the recurrence times of large earthquakes is challenging, even for the simple characteristic\nmodel5-6, mainly due to the unpredictable nucleation conditions7 and the uncertainties in stress states8.\n7. Line 71. What is the geologic estimation of slip rate on DLSF? And how does the geodetic slip rate\ncompare to observations?\nResponse: Thanks for the comment. We explained in the main text:\nThe DLSF has a slip rate of 3.3-4.8 mm/a, responsible for ~40% of the sinistral motion across the fault\nsystem, which is also consistent with the geologic estimation of ~3-4 mm/a23.\n8. Line 77-86. Not sure if earthquake sizes can be inferred from fault maturities. It may be necessary to\ncompare fault systems with similar geologic accumulated slip, age, and tectonic loading condition, and\nearthquake histories to make such calls.\nResponse: We downplay this implication in the main text:\nIf the experimental results of fault roughness in the laboratory (e.g., ref 35) can be upscaled to natural\nfaults, the AZF may tend to produce larger events than the DLSF, yet a conclusive comparison should\nalso consider their accumulated slips, geological ages, tectonic loading conditions, and other relevant\nfactors.\n9. Line 127. Are you simulating the whole Anninghe-Zemuhe fault in dynamic rupture models?\nResponse: We added the following explanation in the main text:\nOn the AZF, we primarily focus on the Anninghe Fault and the northern segment of the Zemuhe Fault,\nas the shear stress on the southern segment of the Zemuhe Fault is currently insufficient to accommodate\nlarge events.\nIf so, how do you justify the planar fault assumption given the change in strike.\nResponse: We have performed additional simulations with curved faults and presented the results in the\nnewly added Figure 4. Please refer to the response to comment 1.\n10. Line 143-144. It would be nice to show rupture extents from the models and compare them to the\ncoupling ratio distribution.\nResponse: We added new Figure S13 to compare the predicted rupture extent with the coupling ratio\ndistribution.\nAre there any multi-segment ruptures comparable to those paleoseismic M 7.5 events?\nResponse: We find all simulated ruptures are confined within each segment in the current stress state\n(Fig. 4C). We expect that ruptures could bridge multiple segments if the AZF remains quiet for hundreds\nof years, which warrants further investigation in the context of earthquake cycle simulation.\nAre there even larger events allowed by the models?\nResponse: We believe it is possible if the AZF remains quiet for a sufficiently long period.\n11. Line 188: the GPS records are critical to the inversion. It would be helpful to talk about, even\nbriefly, uncertainties and qualities of GPS records. For newly deployed stations, it would be helpful to\nhave a table or data file to show the station locations, EW and NS rates and uncertainties for the new data.\nResponse: Many thanks for the comment. We added a new table and more details about the GPS data in\nthe revised manuscript:\nTo better constrain the coupling ratio distribution along the AZF, we built 24 new near-fault GNSS\nstations s in 2018 and surveyed them three times in 2019, 2020 and 2022 (Fig. 1A). In order to reduce\nseasonal effects, field surveys were performed at the similar time of the year (July). The occupation of\neach survey was larger than 72 hours to ensure the stability of the observation. We solve the velocities\non these new GNSS stations using the strategy shown in Diao et al45, and combine them with the existing\nGNSS velocities46 for the kinematic inversion following the method of Wang and Shen46. The solved\nGNSS velocities and associated uncertainties on these new stations are shown in Table S1.\nTable S1. Velocity of 24 newly-developed GNSS stations under Eurasian plate defined by Wang and\nShen46. Ve and Vn are the East and north velocity components under Eurasian plate, dVe and dVn are\nthe corresponding velocity uncertainties (standard deviation).\nStation Lon. (o) Lat.(o) Ve (mm/a) Vn(mm/a) dVe (mm/a) dVn(mm/a)\nZD19 102.270 29.016 5.892 -4.854 2.9 2.6\nZD20 102.313 28.997 15.494 -10.234 1.4 1.2\nZD21 102.226 28.928 14.052 -10.619 1.5 1.4\nZD22 102.310 28.914 12.394 -8.856 0.9 0.9\nZDC2 102.167 28.618 11.201 -11.039 1.0 0.9\nZD24 102.288 28.627 10.923 -8.732 0.9 0.9\nZD25 102.230 28.508 13.139 -9.358 1.0 0.9\nZD26 102.234 28.433 12.861 -10.837 1.4 1.3\nZD27 102.135 28.422 12.315 -12.311 1.0 0.9\nZD28 102.282 28.098 11.805 -11.152 0.8 0.7\nZD29 102.148 28.077 10.856 -9.851 0.8 0.7\nZD30 102.126 27.994 9.486 -10.971 0.8 0.7\nZD31 102.272 27.967 10.302 -10.307 2.6 2.4\nZD32 102.387 27.729 9.408 -10.087 0.8 0.7\nZD33 102.326 27.680 8.992 -10.289 0.9 0.9\nZD35 102.510 27.510 9.173 -10.804 0.9 0.8\nZD36 102.473 27.461 11.416 -9.590 1.0 0.9\nZD37 102.467 27.609 12.905 -8.048 1.1 1.0\nZD38 102.643 27.235 8.938 -10.169 0.8 0.7\nZD39 102.572 27.198 11.883 -11.204 0.8 0.8\nZD41 102.799 26.930 5.699 -13.808 0.9 0.8\nZD42 102.922 26.770 8.303 -9.536 0.8 0.8\nZD43 103.039 26.803 5.761 -4.768 0.8 0.7\nZD44 102.962 26.631 8.963 -8.285 0.9 0.8\nAre there any vertical components?\nWe did not use the vertical components in the model for two reasons: (1) the studied region is dominated\nby horizontal deformation with faults dominated by sinistral strike-slip character; (2) The uncertainties\nof vertical components are larger than those of horizontal components.\n12. For Figure S5, is this block model 2-D or 3-D? Block models are known to suffer trade-offs between\nfault widths and inverted slip rate. Is there such an issue in this study?\nResponse: Our block model is 3D. The trade-offs between fault widths and inverted slip rate occur when\nthe coverage of GNSS stations is sparse around the fault. However, in this study, we have dense near-\nfield GNSS stations and well-covered middle- to far-field GNSS stations, therefore the trade-off is\ninsignificant.\n13. Figure S5, would be nice to have fault names on the maps.\nResponse: We have added fault names on the maps in Figure S7 (previous Figure S5).\n14. Line 200. While Model 1 that ignores the XJHF generates a slip rate on AZF close to geologic\nobservations, is there any estimation of geologic slip rate on the XJHF? If so, how does it compare to\nAZF? Given the close proximity of the two faults, can XJHF be safely ignored?\nResponse: Many thanks for the valuable comment. Geologic and geodetic investigations suggested\nthat the LJ-XJH fault is dominant by sinistral slip with rates of 1 mm/a and 3 mm/a on its northeastern\nand middle-southwestern segments, respectively (ref58,59). Moreover, besides the LJ-XJH fault, other\nsecondary faults, such as the Litang fault lies within the Chuandian block. To check the effect of these\nsecondary faults, we evaluate the effect caused by them and address the issue in the main text:\nMoreover, to test the effect of secondary faults within the Chuandian block (Lijiang-Xiaojinhe fault\nand Litang fault), we forward calculate the crustal deformation due to the motion of these secondary\nfaults using the “buried-dislocation” model54 and associated fault slip rate58,59. As shown in Fig. S9, such\neffect is ~ 0.5 mm/a near the AZF. After removing this effect from the GNSS velocities, we re-run Model\n1 to get the fault coupling distribution (Fig. S9). Not surprisingly, we find that the coupling model\nremains stable and the effect of the secondary fault is negligible.\n15. Line 201. It would be helpful to briefly introduce the geologic slip rate estimations and\nuncertainties for the faults involved.\nResponse: Thanks for the comments. We list the geologic slip rate and related uncertainty of involved\nfault in Table S2, and revised this sentence in the main text:\nModeling results indicate that Model 2 leads to a slip rate of 9.2-5.7mm/a along the AZF (Model 2 in\nFig. S7) that is higher than the geological result of 6.5±1.0 mm/a22. On the contrary, Model 1 that\nneglects the XJHF, generates a slip rate of AZF (6.1-4.9 mm/a) that agrees better with the geologic\ninvestigations.\n16. Line 203. What is the model setup of the viscoelastic model in terms of model domain, fault\ngeometry and historical earthquake slip implementation?\nResponse: Thanks for the comment. We add more detailed description about postseismic viscoelastic\nrelaxation model in the main text:\nFor this purpose, we use the PSGRN/PSCMP code47 with a layered infinite half-space model that\nconsists of an elastic upper crust and two underlying viscoelastic layers representing the lower crust and\nupper mantle (Fig. S5). The elastic and viscoelastic parameters in the model are from ref48 and refs45,49,\nrespectively. The postseismic viscoelastic relaxation effect is driven by coseismic rupture of historical\nearthquakes, of which the parameters are inferred from refs24,50-52 and listed in Table S2. As shown in\nFig. S6, this effect remains less than 0.5 mm/a despite different viscosities being used.\nWhat about the availability of code PSGRN/PSCMP and the viscoelastic model?\nResponse: The PSGRN/PSCMP code is open source that is freely available. We address this point in the\n“Code availability”:\nThe PSGRN/PSCMP code for postseismic deformation modeling and viscoelastic Green’s function\ncalculation, as well as the SDM code for the slip deficit inversions, are available at the public sever of\nGerman Research Centre for Geosciences (GFZ): ftp://ftp.gfz-potsdam.de/pub/home/turk/wang/. The\nopen-source software SPECFEM3D used in the 3D fully dynamic rupture simulations are available from\nthe Computational Infrastructure: https://github.com/SPECFEM/specfem3d. The open-source software\nZMAP used for b value calculations is available at: https://github.com/zmap/zmap.\n17. Line 204. Do you mean Figure S6?\nResponse: Yes, revised.\n18. Line 212. Is s backslip or backslip rate?\nResponse: It is a typo. We have changed “backslip” to “backslip rate” in the revised manuscript.\n19. Line 214. Do you have the trade-off plot between alpha and data misfit?\nResponse: We added it in Fig. S8:\nFigure S8. Trade-off curve with data misfit plotted as a function of normalized model roughness.\n20. Line 215. Would be helpful to lay out equations used to calculate fault coupling ratio from\ninverted slip rates. Also, would be helpful to show briefly how the uncertainties are evaluated.\nResponse: Thanks for the comments. We add more content to introduce the calculation of fault coupling\nratio and related uncertainty estimation in the main text:\nAfter obtaining the back-slip rate, we obtain the fault coupling ratio (r) by 𝑟 = 𝑠𝑏\nwhere s is the fault\n𝑠\nslip rate. We estimate the uncertainty of the inverted slip deficit rate and coupling ratio using the statistic\nmethod47. The details are as follows: we first generate a set of Gaussian noise signals using the\nobservation uncertainty as the standard deviation; Then, we add them to the observation data to get\nseveral sets of disturbed synthetic data; For each of the disturbed data sets, we repeat the inversion\nprocess to get the coupling ratio on each fault patch; These steps were run 500 times in order to estimate\nthe uncertainty of the parameters statistically.\n21. Line 216-224. Seems this paragraph that discusses viscoelastic effect could be merged with\nLine 202-205?\nResponse: Done. See merged contents in the main text:\n𝑉 is accounted for in our model by simulating the postseismic viscoelastic effect of historical\n𝑃\nearthquakes since 1480. For this purpose, we use the PSGRN/PSCMP code47 with a layered infinite half-\nspace model that consists of an elastic upper crust and two underlying viscoelastic layers representing\nthe lower crust and upper mantle (Fig. S5). The elastic and viscoelastic parameters in the model are from\nref 48 and refs45, 49, respectively. The postseismic viscoelastic relaxation effect is driven by coseismic\nrupture of historical earthquakes, of which the parameters are inferred from refs24,50-52 and listed in Table\nS2. As shown in Fig. S6, this effect remains less than 0.5 mm/a despite different viscosities being used.\nBeside the viscoelastic relaxation effect due to stress loading of historical events (𝑉 in Eq. (1)), we also\n𝑃\nconsider the viscoelastic effect caused by the interseismic stress accumulation due to the fault coupling.\nBased on the theoretical derivation shown in ref53, we consider the interseismic viscoelastic effect by\nconstructing the viscoelastic Green’s function of 𝐺(𝑥,𝑡)= [𝐸(𝑥)+𝑉𝑖𝑠(𝑥,𝑡)𝐻(𝑡)] for slip deficit\ninversions (Eq. (2)), where 𝐻(𝑡) is the Heaviside function, 𝐸(𝑥) the elastic displacement, and\n𝑉𝑖𝑠(𝑥,𝑡) the time-dependent displacement induced by the interseismic viscoelastic relaxation.\nConsidering the long earthquake recurrence period on the AZF, we calculate 𝐺(𝑥,𝑡) by simulating a\nfully relaxed viscoelastic effect (𝑡 → ∞) using the layered viscoelastic model (Fig. S5).\n22. Line 257. It seems a planar fault geometry is used. Is the whole Anninghe-Zemuhe fault system\nsimulated in dynamic ruptures? The fault strike variations are not small to be ignored. What is the total\nfaut length?\nResponse: Please refer to our response to comment 1.\n23. Line 262. What is the CFL coefficient and time step size? Would be helpful to have a table of\nkey parameters.\nResponse: We added a new table for these key parameters and revised the text in Methods:\nThe time step is set as 0.007s, meeting the Courant-Friedrichs-Lewy stability condition, a criterion\nwidely used to ensure the stability of the numerical simulation.\nTable S3. Key model parameters used for dynamic simulation models in this study.\nPower-law slip weakening Exponential slip weakening\nComputational domain 200 km *100 km*70 km 200 km *100 km*70 km\nTime step 0.007s 0.007s\nAverage grid size 250 m 250 m\n𝜏 −𝜏 28 MPa 8 MPa\n𝑠 𝑑\nd 0.01 m 0.2 m\nc\np 0.15 -\n24. Line 265. Same to the comment to Line 23.\nResponse: Please refer to our response to comment 4.\n25. Data Availability. Is there any public data repository to store and share the data?\nCode availability. Nature (and other journals in general) recommends DOI-minting repository to store\nand share code. What about the code availability of PSGRN/PSCMP? Is Computational Infrastructure\na data repository?\nResponse: We add required information in the 'Code availability':\nThe PSGRN/PSCMP code for postseismic deformation modeling and viscoelastic Green’s function\ncalculation, as well as the SDM code for the slip deficit inversions, are available at the public sever of\nGerman Research Centre for Geosciences (GFZ): ftp://ftp.gfz-potsdam.de/pub/home/turk/wang/. The\nopen-source software SPECFEM3D used in dynamic rupture simulations is available from the\nComputational Infrastructure: https://github.com/SPECFEM/specfem3d. The open-source software\nZMAP used for b value calculation is available at: https://github.com/zmap/zmap.\nand in 'Data availability':\nThe velocity of newly-developed GNSS stations is listed in Table S1, and the other GNSS velocities can\nbe found in the supplementary files of ref46. The refined seismic catalog used for b value calculation is\nshown in Supplementary Data 1.\nREVIEWERS' COMMENTS\nReviewer #2 (Remarks to the Author):\nI appreciate the authors for their diligent efforts in addressing my comments. The\nmanuscript is much improved. I endorse its publication now. Thanks!\nOne last comment: I still feel there is room for improvement related to my previous\ncomment \"In Figure 2a, ruptures that connect two high coupling segments (1480 M7.5, 1536 M7.5, 1850 M7.5)\nare mixed with ruptures limited within single high coupling segment (1732 M6.8, 1952 M6.8, and\npotentially 1489 M6.8). There are also interesting dynamics that the 1489 M 6.8 was followed shortly\nbe the 1536 M7.5 rupture, and that the 1732 M 6.8 rupture was followed by the 1850 M 7.5\nrupture.\nDynamic rupture model is an adequate tool to study multi-segment ruptures.\"\nThese potential multi-segment ruptures are relevant to seismic hazard evaluation. I look\nforward to seeing progress from your method to address the potentials of multi-segment\nruptures, and the complex dynamics in terms of both their rupture extents and recurrences.\nFault geometric complexity affects rupture extents through heterogeneous stress conditions\nthat are inherited from past earthquakes and developed over cycles. In Figure 4C, it appears\nthat ruptures are well confined within bend limits. So L184 may not be true. Anyhow, they\ndon't contradict to your method, because the heterogeneous stress states are well\ncharacterized through your inversion with GPS constraints. It is nice to see a planar fault\nperforming so well with your method.\nResponse to reviewers' comments: Black Italic = Comment, Blue = Response,\nRed = Revision\n=============================================================\nReviewer #2 (Remarks to the Author):\nI appreciate the authors for their diligent efforts in addressing my comments. The manuscript\nis much improved. I endorse its publication now. Thanks!\nResponse: Thank you for accepting our revisions.\nOne last comment: I still feel there is room for improvement related to my previous comment\n\"In Figure 2a, ruptures that connect two high coupling segments (1480 M7.5, 1536 M7.5, 1850\nM7.5) are mixed with ruptures limited within single high coupling segment (1732 M6.8, 1952\nM6.8, and potentially 1489 M6.8). There are also interesting dynamics that the 1489 M 6.8 was\nfollowed shortly be the 1536 M7.5 rupture, and that the 1732 M 6.8 rupture was followed by\nthe 1850 M 7.5 rupture. Dynamic rupture model is an adequate tool to study multi-segment\nruptures.\"\nThese potential multi-segment ruptures are relevant to seismic hazard evaluation. I look\nforward to seeing progress from your method to address the potentials of multi-segment\nruptures, and the complex dynamics in terms of both their rupture extents and recurrences.\nFault geometric complexity affects rupture extents through heterogeneous stress\nconditions that are inherited from past earthquakes and developed over cycles. In Figure 4C,\nit appears that ruptures are well confined within bend limits. So L184 may not be true. Anyhow,\nthey don't contradict to your method, because the heterogeneous stress states are well\ncharacterized through your inversion with GPS constraints. It is nice to see a planar fault\nperforming so well with your method.\nResponse: Thank you for this very interesting comment. Due to the lack of slip models for\nhistorical events, our current model cannot investigate potential connections between those\nsmall and large events. The topic of potential multi-segment ruptures will be the focus of our\nfuture work, which will combine the rupture theory, long-term earthquake cycle simulation, and\nseismological observations. To acknowledge this discussion, we have added one sentence in\nthe main text:\nBut its neighboring asperities can still accommodate Mw > 7 events and their maximum\nmagnitudes slightly increase due to the loading from the 1952 rupture. Due to the lack of direct\nconstraints on the 1952 Mw 6.8 earthquake in the current model, the possibility of a future Mw\n7.5 event bridging AFZ-3 and AFZ-4, similar to the 1536 M 7.5 event, cannot be excluded and\nthus is worthy of further investigation.\nIn addition, we agree that the constrained heterogeneous stresses may stem from fault geometric\ncomplexity. We revised L184 to make it clearer:\nBased on the current stress state of the AZF constrained in this study, we find fault bending is\nnot the primary factor controlling earthquake sizes, ..."}]

Decision letter