A systematic scheme to develop dynamic earthquake rupture scenarios: a case study on the Wenchuan–Maoxian Fault in the Longmen Shan, China, thrust belt

The 2008 Wenchuan Mw 7.9 mainshock caused catastrophic destruction to cities along the northwestern margin of the Sichuan Basin. This earthquake did not activate the Wenchuan–Maoxian Fault (WMF) on the hinterland side and the conjugate buried Lixian Fault (LXF), but they could experience large earthquakes in the future. We propose a systematic scheme to develop scenario earthquakes for active fault systems with insufficient constrain of 3D fault geometries. We first performed stress tensor inversion to constrain the regional stress field. Then, we developed a new method to constrain fault geometries by inverting long-term slip rates under the given regional stress and applied it to the WMF. We conducted a set of 3D dynamic earthquake rupture simulations on the WMF and LXF to assess the scenarios of earthquake rupture processes. Several fault nucleation points, friction coefficients, and initial stress states are assessed, the general rupture patterns for these earthquake scenarios are evaluated, and finally, we find the scenarios that could fall into three groups. Depending on initial conditions, the dynamic rupture may start in the LXF, leading to magnitude-7.0 earthquakes, or start in the WMF, then cascade through the LXF, leading to mag-nitude-7.5 earthquakes, or both start and arrest in the WMF, leading to around magnitude-6.5 or -7.0 earthquakes. We find that the rupture starting on the reverse oblique-slip jumps to the strike-slip fault, but the reverse process is impeded.


Introduction
The 2008 Mw7.9 Wenchuan earthquake caused devastating destruction to cities and counties along the active faults of the Longmen Shan thrust belt.Three parallel NW-dipping fault zones from the hinterland to the foreland [the Wenchuan-Maoxian Fault (WMF), the Beichuan Fault (BCF), and the Pengguan Fault (PGF)] constitute the crustal structure of the Longmen Shan Fault Zone (LMSFZ) (Fig. 1).The mainshock and aftershocks only ruptured the BCF and PGF (Xu et al. 2009;Feng et al. 2017;Liu-Zeng et al. 2009).Several counties, including Yingxiu, Beichuan, and Nanba, suffered a disaster in this event.One reason is that the fault almost passes through the center of these counties (Fig. 1).Similarly, the WMF also passes through the urban areas of Wenchuan, Maoxian, and several towns.Nearly 200,000 inhabitants live in the area along this fault.Should an earthquake in the future occur on the WMF, it would cause a significant hazard to densely populated towns.Assessing earthquake rupture scenarios in the future is essential to increase society's preparedness.
Dynamic rupture simulations intrinsically simulate physically self-consistent features of earthquake behavior and reproduce retrospectively large earthquake events (Aochi and Fukuyama 2002;Ando et al. 2017;Ando and Kaneko 2018;Ulrich et al. 2019).Since the dynamic rupture simulation is a forward problem to constrain the model based on the pre-seismic observations, efforts are made to use it to explore future scenario earthquakes under pre-seismic estimation of stress conditions and fault geometries (Oglesby and Mai 2012;Zhang et al. 2017;Aochi and Ulrich 2015;Ramos et al. 2021;Harris et al. 2021;Yao and Yang 2022).
Three-dimensional fault geometry is a dominant factor in earthquake rupture dynamics.Irregular fault geometry, such as bends, branches, intersections, and stepovers, have been evidenced to cause rupture initiation and termination (Oglesby 2005;Elliott et al. 2015;Yang et al. 2013;Ando and Kaneko 2018).Several studies (Shi et al. 1998;Oglesby et al. 2000;Ma and Beroza 2008;Tang et al. 2021a) have shown that the dip angle, which leads to the asymmetric geometry of dip-slip or oblique faults intersecting the free surface, can affect the dynamics of earthquake rupture and ground motions considerably.In earthquake simulations, fault geometry is typically obtained from detailed relocation of aftershocks (Wu et al. 2009;Zhao et al. 2011;Yin et al. 2018), inversions using interferometric synthetic aperture radar (InSAR) measurements (Shen et al. 2009;Wan et al. 2016;Fukahata and Hashimoto 2016) and seismic reflection data (Hubbard et al. 2010;Jia et al. 2010;Lu et al. 2014).These methods are only available to investigate earthquakes that have already occurred.Seismic reflection and borehole surveys can estimate pre-seismic estimates of fault geometry, but they are costly in operations, and the data are not always available.Additionally, gravity and magnetic (Tian et al. 2017), seismic tomography (Tong et al. 2019), and magnetotelluric sounding (Zhao et al. 2012) cannot provide enough spatial resolution to describe the fault geometry, especially the dip angles.More straightforward and more effective methods have been awaited.
We propose a new method to infer the fault geometry by inverting the long-term slip rate based on the Wallace-Bott hypothesis (Wallace 1951;Bott 1959) to refine the fault geometry based on the existing inference of the WMF (Hubbard et al. 2010) and the surface fault trace (Xie et al. 2011).This hypothesis recalls that two bodies in contact slip toward the interfacial tangential traction and asserts that faults slip in the same manner (Fig. 2).This is also the base of the widely used stress tensor inversion method (Hardebeck and Michael 2006).
Moreover, this hypothesis is confirmed based on a large natural earthquake's seismically inverted slip distributions (Matsumoto et al. 2018).Recent detailed analyses of InSAR observations demonstrate that the slip direction at the ground surface is consistent with those of the overall characteristics at seismogenic depths (Fukahata and Hashimoto 2016).Therefore, in an ideal case, the fault plane and the principal stresses should uniquely determine the sense and direction of the long-term slip rate.Since geological and geomorphological observations can determine fault strikes, principal stress axes, and longterm slip rates, the Wallace-Bott hypothesis suggests that dip angles can be inferred based on these observations.In fact, this idea is a natural and quantitative extension of our commonly used experience to distinguish the near-vertical strike-slip faults from more tilted dip-slip faults.An advantage of this method is that we can provide observation-based and systematic inferences of the fault geometry even if more direct data such as seismic reflections survey and hypocenter locations are absent.
In this study, we first calculate the distribution of heterogeneous stress fields near the WMF based on the focal mechanism solutions (FMS) of 2008 Wenchuan earthquake aftershocks (M ≥ 3.5).The heterogeneous on-fault initial stresses are resolved from the tectonic stress tensors.Then, we use the long-term slip rate together with the inferred stress to refine the WMF structure based on the Wallace-Bott hypothesis.We then implemented dynamic rupture simulations to investigate rupture propagations on the WMF and surrounding faults to assess potential rupture scenarios.(Xu et al. 2009;Liu-Zeng et al. 2009).The blue solid lines represent the WMF and LXF, both of which did not slip in the mainshock but are very active in aftershocks.The green solid lines represent three active faults parallel to the LXF.Accurate fault trace of the WMF is derived from field investigation of Xie et al (2011), and the location of the LXF trace is deduced from the distribution of pure strike-slip aftershocks (Li et al. 2019).Historical earthquake catalog denoted by filled circles refers to Wang et al (2021) and China Earthquake Networks Center (http:// data.earth quake.cn).Green hollow circles denote the historic strong earthquakes (M > 5) before 2008 Wenchuan event.The purple and gray coupled arrows depict the maximum and minimum horizontal compressional stress fields (Luna and Hetland 2013).The third primary stress (not shown) is vertical.The purple dashed regions, respectively, show the positions of Xuelongbao (north) and Pengguan (south) Massif (Shen et al. 2019)

Seismic and tectonic condition of the Wenchuan-Maoxian fault and Lixian fault
The WMF has prominent fault orientation and geometry similarities to the BCF and PGF that caused the 2008 mainshock.The three faults are subparallel and horizontally spaced within 20 km (Fig. 1) and might root into the same main detachment with a depth of approximately 20 km.The main fault shapes display a classic ramp-flat geometry (Hubbard et al. 2010;Jia et al. 2010;Li et al. 2010).Tang et al. 's (2021b) simulation revealed that the Wenchuan mainshock changed the stress on the WMF, generating a positive Coulomb-failure stress change ( CFS ) in the center portion of the WMF, suggesting that the 2008 mainshock might have clock-advanced the next earthquake on the WMF to some extent.The WMF has not experienced a strong earthquake for more than 300 years since the last M-6.5 earthquake in 1657 (National Earthquake Data Center, http:// data.earth quake.cn), and the current earthquake interval might provide sufficient accumulated strain for a moderate earthquake (5.5 < M < 6.5) (Additional file 1: Fig. S1).
Field studies show that the WMF primarily exhibited dextral strike-slip faulting with a comparable reverse component (Tang 1991;Rongjun et al. 2007;Shen et al. 2019;Tian et al. 2013;Tan et al. 2017) since the Late Cenozoic, consistent with the coseismic slip sense and the long-term slip rate of BCF and PGF (Rongjun et al. 2007;Densmore et al. 2007;Xu et al. 2009;Liu-Zeng et al. 2009;Feng et al. 2017).Based on fission track and (U-Th)/He dating, Shen et al. (2019) obtained a 0.6 mm/ yr long-term slip rate of WMF, which is as high as those along the frontal BCF (~ 0.54 mm/yr before the Wenchuan earthquake (Rongjun et al. 2007) and 0.88-0.91mm/yr after the Wenchuan event (Ran et al. 2013)].
In addition, according to observational seismic records before the Wenchuan earthquake (Additional file 1: Fig. S1), four earthquakes occurred near the WMF out of 11 moderate earthquakes (5 < M < 6) in this area.Note that 12 Ms > 5.6 aftershocks occurred within 1 year after the 2008 Wenchuan earthquake, of which one Ms 6.1 event occurred on the southernmost end of the WMF and one Ms 6.0 event in the Lixian Fault (LXF) (Zhao et al. 2011).
The buried LXF is geometrically close to the WMF and bridges between the WMF and the BCF by its conjugate geometry (Fig. 1), which may result in cascading ruptures to form larger earthquakes.Geological studies are scant, and the activity is yet to be clarified for this fault due to the heavy vegetation cover and rugged terrain.The maximum altitude difference exceeds 3 km (Tan et al. 2017).The LXF appears to extend to the Miyaluo Fault, a welldefined geological structure (Fig. 1) displaying a sinistral sense of slip with unknown slip rates and thought to be active at least during the Pleistocene based on the fission track dating (Yang and Zhang 2010).Several earthquakes larger than M-7 were recorded on the Fubianhe and Songpinggou Faults, which are subparallel to the Miyaluo Fault with an NW strike orientation (Fig. 1).In the following section of the stress inversion analysis, the LXF exhibits the most favorable orientation in the regional stress field.Therefore, the LXF might display a comparatively high seismic potential, although none of the great earthquakes has been recorded to have occurred on it recently.

Regional stress inversion and modeling
The initial stress conditions are fundamental in controlling dynamic earthquake rupture processes (Harris 2004;Ramos et al. 2022).We inferred the present regional tectonic stress field by inverting earthquake Fig. 3 Along strike distribution of three-component principal stress ( σ H ,σ h for the maximum and minimum horizontal principal stresses, and σ v for vertical principal stresses).The lower right figure shows the angle between the east and the maximum horizontal stress axis (negative denotes the clockwise), and black box denotes the projection of LXF FMSs (from Li et al. 2019) of 2008 Wenchuan earthquake aftershocks for computationally simulating earthquake scenarios in this region.We adopted a damped linear stress inversion method developed by Hardebeck and Michael (2006) to determine the stress field using FMSs.This method enables us to infer the directions of the principal stresses σ r (r = 1, 2, and 3 for the maximum, intermediate, and minimum principal stresses, respectively) and the ratio between the principal stresses σ r , defined as ∅ = (σ 2 − σ 3 )/(σ 1 − σ 2 ) .The details of the method and results are described in Additional file 1: Text S1.
However, in our dynamic rupture simulation, these observationally constrained parameters are insufficient to determine the absolute values of the principal stresses.Two unknowns are needed and obtained by imposing other prior information.For simplification, we first assume that one of the principal stress axes is vertical based on the observation near the WMF, as by Anderson (1905).Then, we constrain σ H /σ v to 1.6, an approximate value from the in situ stress measurements along the BCF after the Wenchuan earthquake (Qin et al. 2015) and previous simulation works (Duan 2010;Zhang et al. 2019;Tang et al. 2021b), where σ H and σ v denote the magni- tudes of the maximum horizontal and vertical principal stresses.The σ H along the depth is constrained by assum- ing that the average stress drop satisfies the stress drop derived from empirical relations (Kanamori and Anderson 1975).We assume typical overburden stress and that the pore pressure increases with the depth, and the latter gradient gradually approaches the former (Rice 1992), so that the stress is tapered to the ambient value at a specific depth (here, we set 5 km) (Rice 1993).Figure 3 shows the resulting stress field model of the WMF.The along-strike rotation of the principal horizontal stresses reflects our inversion result.Since there are no independent constraints on the specific depth of the WMF, we discuss its influence in the section "Effect of saturation depth of stress".

Fault geometry construction
Our dynamic model must match the observed oblique dextral slip WMF and the almost pure sinistral strikeslip LXF in the long-term slip rate.We set the LXF as vertical, with a 0.5 km buried depth and a 15 km basal depth inferred from the average seismogenic depth (Shen et al. 2009;Wan et al. 2016;Ramirez-Guzman and Hartzell 2020).In the simulation, we model only the southern portion of the LXF, and do not consider its possible extent further to the north, which could be hypothetically a link to the Miyaluo Fault (Yang and Zhang 2010).
The aftershocks of the M-7.9 Wenchuan earthquake and the 2013 M-7.0 Lushan earthquake delineated a prominent seismic gap approximately 50 km long near the town of Dayi (the blue box in Fig. 1).Several studies suggest that this gap can correspond to a possible unruptured fault segment (Wang et al. 2015(Wang et al. , 2018;;Liu et al. 2018), and, in this study, we assume that the south part of the WMF terminates near this gap.
Although we could mimic this fault system with several planar rectangular fault segments, we developed a nonplanar fault model by considering existing data and geological inferences.We developed the WMF geometry model by considering (1) geological investigations and seismic reflection data (Xu et al. 2009;Hubbard et al. 2010;Jia et al. 2010), providing an approximate geometry of the WMF at a depth: a fault subparallel to the BCF, with a larger dip angle and listric geometry.(2) The surface fault trace from field investigations (Xie et al. 2011) provides along-strike constraints on the surface fault structure.(3) The early aftershock locations (Yin et al. 2018) help constrain the south section of the WMF geometry.Based on these constraints, we refine the geometry by developing a new method by inverting the long-term slip rate and the regional stress field based on the above consistency between the slip direction and the tangential traction (Wallace-Bott hypothesis) (Fig. 2).The direction of the on-fault tangential traction is fully determined from the stress tensor inversion analysis, including the orientations of the three principal stress axes and the ratio between the magnitudes of the principal stresses (stress ratio) (Michael 1984).The validity of this new method is observationally supported by Matsumoto et al (2018); in the case of the 2016 M7.0 Kumamoto, Japan, earthquake, they demonstrated the good coincidence of the directions between the coseismic slip The readers should not confound our inferring dip angles with the internal or optimal frictional angles determined by the frictional coefficient.The dip angle inferred in our method is independent of the friction coefficient.In other words, our method does not hypothesize that faults occur in the most optimal orientation and has the flexibility to apply to misoriented faults under the recent tectonic stress conditions.Our approach can be valid if the regional stress fields are stable over time to have the observed long-term slip rates.While the temporal stability of tectonic stress regime (or the absolute stress level, not the stress change) considered here is an issue that needs more intensive investigations, at least, it appears to be insensitive to individual occurrence of large earthquakes (i.e., Imanishi et al. 2012).

Geometry inversion method
A proper choice of geometry models should be made depending on purposes by considering the available data and previous inferences (Xu et al. 2009;Hubbard et al. 2010;Jia et al. 2010).A model of listric fault geometry is considered here by following the above models.The general dip angles are allowed to vary continuously along the fault strike, and the surface fault traces are fixed, as given by the data.We use one inversion parameter η to describe the fault dip in the local area by assuming the listric fault geometry, as described by (Wan et al. 2016) where z is the depth of the fault plane, and y is the horizontal distance of the fault plane from its surface trace.The fault surface is described by the B-spline varying continuously along the strike, with its downdip curvature dictated by the parameter η (Fig. 4).h 0 is the fault depth to invert and is assumed to be the same for all fault segments.n B-spline curves are uniformly distributed along the fault strike with equal horizontal intervals X to The direction of the long-term slip rate is approximately collinear with the traction on the fault plane, which is equivalent to solving the problem of minimizing the difference between the two ratios.We solve the inverse problem in the sense of Tikhonov (1963), taking a regularized solution to be a model minimizing an objective function, , defined by ( 3) where T 1 and T 2 denote the averaged traction along the strike and dip within each spherical computational region with radius R, respectively, represented at the measured points.The averaged traction introduced here can be regarded as a comprehensive effect within the computational region R and depends on the regional stress and fault geometry (Fig. 4).The V 1 and V 2 denote the hori- zontal and vertical long-term slip rates at the measured points.We choose the damping matrix L to be a simple, first-difference operator (make sense for only η not h ).
The regularization parameter is a positive number that balances the effects of data fitting and model smoothness, controlling the smoothness of model space.However, the observed slip-rate ratios are not always obtained for each measured point, where either horizontal or vertical long-term slip rates can be given.We prefer to use an approximate form where α is a scaled value manually set to meet the best initial misfit.We construct the forward model to connect the fault geometry parameters to weighted average tractions T 1 (m) and T 2 (m) near the surface (Fig. 4).We used the conjugate gradient method to obtain the optimal solution of m .The forward calculation and inversion pro- cedures were further elaborated in detail in Additional file: 1 Appendix A of the supplementary material.

Application to the WMF
In this study, the model parameter m include ten down- dip curvature parameters η uniformly distributed along the fault strike, plus h 0 for all curves.Under the condi- tion that the values of η can accurately describe the sur- face trace and the nonplanar fault surface, we expect as few parameters as possible to avoid the overfitting problem due to limited observed data.We collected the vertical and horizontal components of the long-term slip rate from a few studies (Ma et al. 2005;Rongjun et al. 2007;Ran et al. 2013;Shen et al. 2019) shown in Fig. 5 and implemented five measured points with seven slip rates along the WMF in the inversion.We choose the B-spline interval X = 16 km, radius R = 15 km, and the update strategy of the regularization parameter based on the additional tests (Additional file 1: Text S2).
The long-term slip rate cannot be accurately measured due to the uncertainty of the geological age estimation of sediments (Rongjun et al. 2007;Shen et al. 2019).The inversion uncertainty is estimated using 100 bootstrap resamplings of the entire data set (seven slip rates) to obtain a more reliable result.The data from each measured point obey a Gaussian distribution. (4) For each inversion, we set the same initial model as m = [η 0 , . . ., η 0 , h 0 ] , where η 0 = 5.0 , h 0 = 18 km (Fig. 6a).
A few aftershocks are distributed at the south end of the WMF (Additional file 1: Fig. S2), which is helpful as prior information to constrain one of the inversion parameters.
We can fit the aftershock distributions using Eq. ( 1) with η 0 = 5.0 , h 0 = 18 km.Therefore, the parameter η at the southernmost position (x = 0 km) is enforced for little variation during the inversion.An average of 100 inversions is derived to obtain the final fault model (Fig. 6b), where the dip angles near the surface vary along the fault strike due to the constraint of the observed long-term slip rates (Fig. 7a).The final inversion results are distributed in a narrow strand (Fig. 7b).All objective functions steadily decrease versus the iterations (Fig. 7c) and uncertainty of inverted dip angle is within 10° (Fig. 7d), justifying our method.The predicted average fault dip angles at the intersection of WMF and free surface are approximately 60°, typically larger than that in the southern section of the BCF.The above results are consistent with the previous studies (Xu et al. 2009;Hubbard et al. 2010;Jia et al. 2010).Differences in the dip angle can lead to differences in rupture scenarios, including the strong ground motion and surface offset for the BCF and WMF.
While we believe that this method of fault geometry inversion is an effective physically based and observationally supported tool, this method could involve ambiguities due to the limitations of the available data and model constraints in inversion analyses.Furthermore, the Wallace-Bott hypothesis is not always satisfied, as the slip direction might deviate from the tangential traction orientation by considering several factors, including local variation in fault geometry (Lejri et al. 2015;Lisle 2013).Nevertheless, we can find the utility of our optimized search method when we wish to refine previous models (Hubbard et al. 2010) or decide on a better model based on prior information without more direct geophysical or geological information.

Method
We computationally simulate spontaneous dynamic rupture propagation on the conjugate fault system of the WMF and LXF.The dynamic rupture problems are numerically solved by the Fast Domain Partitioning Boundary Integral Equation Method (FDP-BIEM) (Ando 2016;Ando et al. 2017;Ando and Kaneko 2018), which increases the efficiency of the elastodynamic boundary equation method (BIEM) without degradation in accuracy.We introduce additional 'free surface elements' , on which a free surface boundary condition is satisfied (Ando and Okuyama 2010;Hok and Fukuyama 2011).This method could avoid further complicated theoretical computations for the integration kernels, and enable us to introduce rather complicated geometry of fault systems in a half-space medium.We consider a homogeneous elastic half-space, with Lamé parameters λ = μ = 28 GPa and mass density ρ = 2776 kg/m 3 (v p = 5500 m/s, v s = 3175 m/s).These parameters are based on the firstorder approximation of a three-dimensional velocity structure around the Longmen Shan Fault obtained from seismic tomography (Pei et al. 2010;Wang et al. 2021).The local variation of shallow seismic velocity structures or site effects amplifying ground motion is excluded from this model.Instead, we focus on the source effect.
The boundary condition for the fault surface is considered where the frictional strength is described by the linear slip-weakening friction law (Andrews 1976) (5) , where D c , µ s , and µ d denote the characteristic slip and the static and dynamic frictional coefficients, respectively.τ and σ n represent the shear traction and normal stress, respectively.Based on the results of high-speed friction experiments (Yao et al. 2013), we set the friction coefficients uniformly on the entire fault with μ s = 0.53 and μ d = 0.12.We choose D c = 0.8 m, because its upper bound is estimated to be 1.0-2.0m by Tang et al. (2021b), who applied the approach proposed by Mikumo et al. (2003) to estimate the upper bound of D c of the 2008 Wenchuan earthquake, constrained by the slip rates at fault locations from the kinematic inversion (Zhang et al. 2014) of near-field waveforms.The neighboring parameter space is also investigated and discussed in the section "Discussion", considering the uncertainties of these parameters.We first chose moderately low μ s , μ d , and D c to allow the rupture to propagate to the artificially assumed boundaries without being arrested to explore large earthquake scenarios.As for the numerical model size, we typically employed the boundary elements of approximately 26,824.We control triangular elements to be 500 m on the fault and free surface in the vicinity of the fault, and then increase this size with a maximum rate of 7.5 times toward the outer boundary of the free surface until it reaches the maximum size of around 2000 m, such that higher calculation accuracy can be achieved near the fault plane, and computational expense can be greatly reduced due to the gradual distribution of mesh density (Additional file 1: Fig. S3).
In the simulation, the ruptures are triggered by overstressing a small circle patch with a radius of 3 km , where the initial stress is uniform and slightly above the yield stress.Once a rupture is triggered, its subsequent development is controlled by the elastodynamic equation, the friction law, the initial stresses, as well as the fault geometry.Under our current state of knowledge, predicting a likely hypocenter on the fault system of the WMF and LXF is challenging.Therefore, we considered five nucleation locations at the center, the east and west ends of the WMF, and the north and south ends of the LXF.Considering stress changes caused by Wenchuan mainshock and aftershocks (Additional file 1: Text S3), we set one of the potential epicenters near the Wenchuan County.

Results
Figure 8 shows the total slip distribution on the WMF and LXF and the moment magnitude for each case.The rupture tends to jump from the oblique-slip fault (in the WMF) to the strike-slip fault (in the LXF), but the reverse is more difficult.For instance, a dynamic rupture might start in the central, eastern, or western parts of the WMF, and then cascade to the LXF, leading to three comparable-sized earthquakes with M-7.56, M-7.55, or M-7.57.However, if it starts in the south or north of the LXF, triggering a second rupture on the WMF will be challenging (Fig. 8, video S1-S2), leading to a smaller earthquake with moment magnitudes 7.05 or 7.06.
The total slip on the LXF is smaller than that on the WMF.This might be surprising, because the stress drop on the LXF is significantly larger than that on the WMF (Additional file 1: Fig. S4).Noting that the friction coefficients on the two faults are the same, and the faulting and geometry should cause the prominent difference between the coseismic slips.The WMF is a thrust fault with a larger extent in dip and strike than the LXF, which is a buried and smaller purely strike-slip fault.Under the same stress conditions, larger faults tend to slip more, because healing phases suppressing slip may arrive later, and the free-surface effect is reduced for buried faults.Numerical simulations show that thrust faults have larger coseismic slips than normal and strike-slip faults (Oglesby et al. 2000;Tang et al. 2021a) named the hanging wall effect (Abrahamson and Somerville 1996), because the reflected waves from

Table 1 Models with a different set of friction and nucleation positions in this study
H of σ v *: the saturation depth of the two different depth-dependent stress cases shown in Fig. 9 Arrest*: if Yes, the rupture will be arrested rather than reaching the assumed artificial barrier Dynamic triggering*: if Yes, the rupture on one fault will trigger the secondary rupture on the other fault The first letter of the model's name represents the nucleation location, and E1 ~ C1 have been elaborated in the result section the free surface amplify the motions of the reverse fault near the free surface.Therefore, reverse faults tend to slip more than strike-slip faults, even in relatively lower stress conditions.
Rupture directivity effects caused by different nucleation positions lead to a prominent spatial difference, confirmed by surface displacements and PGV (Additional file 1: Text S4).The site toward which the rupture front propagates will suffer stronger ground motions than the site where the rupture propagates away due to the cumulative effect of the seismic radiation called forward directivity (Somerville et al. 1997).We further evaluate the strong ground motion and surface offset in Additional file 1: Text S4.

Discussion
We explored several large earthquake scenarios.However, the dynamic rupture propagation results depend on several key factors, including the initial stress, fault geometry, frictional parameters, and hypocenters, which might lead to large uncertainty for different earthquake scenarios.For instance, different cases with a linear increase in stress from the surface down to some specified saturation depth might result in different ground motion distributions.We include this uncertainty by changing the specified depth but lowering the maximum stress, so that the same average stress drop can be maintained for the two cases (Fig. 9).Without loss of generality, we model 20 scenarios of dynamic rupture propagations (Table 1) to investigate the effects of neighboring parameter space, including μ d (0.12-0.15) (Yao et al. 2013), D c (0.8-1.2 m) (Tang et al. 2021b), and the stress saturation depth.

Effect of saturation depth of stress
The average stress drops on the WMF and LXF are obtained from observations from the 2008 Wenchuan earthquake, but the distribution of stresses along the depth is poorly constrained.The magnitude of the principal vertical stress depends on the overburden pressure (density and depth), which might also be compensated for with increased pore pressure (Hardebeck and Okada 2018).Thus, we build two depth-dependent shear and normal stress regimes with the same average stress drop �σ w = 2.9 MPa.The result of the first case has been discussed in the result section.However, the second case leads to different earthquake scenarios and distributions of the coseismic slip (Fig. 9).For the second case with deeper saturation depths, simulated ruptures starting in the east or center of the WMF propagate to 20 km beyond the nucleation area and arrest at the patch where the normal stress is stronger (Additional file 1: Fig. S4).Furthermore, the rupture from the west part propagates over the WMF and cascades to the LXF but leads to a lower magnitude (Table 1) and coseismic slips than that of W1 (Fig. 8).However, unlike the scenarios on the oblique fault (the WMF), the coseismic slip and magnitudes on the sinistral fault (the LXF) are slightly larger than that in case 1.These phenomena suggest that depthdependent stress patterns might affect dip-slip faults more than strike-slip faults, and a larger stress gradient near the free surface of an oblique-slip fault is more likely to trigger a stronger coseismic slip.

Patterns of ruptures in the WMF and LXF fault systems
The earthquake scenarios from all simulations in Table 1 can be summarized as six patterns (Fig. 10a).The 20 cases evaluated do not cover the entire range of possibilities, and the assigned probabilities remain arbitrary, because the uncertainty of the geometry and initial stress is not fully considered.However, the proposed scenarios in this study are understandable from a mechanics standpoint; therefore, they appear possible.The spatial heterogeneity introduced in this study originates from fault geometry and inverted regional stress, which principally controls the macroscopic patterns.For instance, cascading ruptures jump from oblique to strike-slip faults, but the reverse process is less likely.The ruptures are arrested on the WMF, whereas they easily run through the entire LXF (Fig. 10a).Furthermore, a rupture propagating Figure 10b represents a statistical analysis of the moment magnitude based on the deterministic simulation results.The earthquake scenarios fall into three groups: magnitudes of approximately 7.5, 7.0, or near 6.6.The first group includes earthquakes that start on the WMF and can propagate over all or most parts of the two faults.Earthquakes with magnitudes of approximately 7.0 typically come from the LXF or local ruptures on the WMF.The third group includes earthquakes that start on the WMF but quickly stop near the nucleation patch due to the unfavorable fault orientation because of the friction parameters and on-fault stress.
Incorporating all earthquake scenarios to further assess the potential seismic motion, we calculate the maximum surface displacements shown in Fig. 11.Static surface deformation is crucial for near-source hazard analyses.Figure 11 shows the two prominent peaks at 40 km and 100 km along-strike distances.Wenchuan town is near the edge of the moderate slip, and Maoxian town is in the extent of the larger slip.The PGV distributions (Additional file 1: Text S4) indicate moderately higher values of PGV around Maoxian than Wenchuan when a rupture nucleates at the center or eastern end of the WMF.The velocity structures, including site conditions, are not included in these simulations; therefore, the results for ground shaking from an actual earthquake, even if it were the same source as shown in these simulations, would differ significantly due to the structural effects.
Furthermore, the values of μ d = 0.12-0.15are slightly lower than those used in Wenchuan earthquake simulations (μ d = 0.18) (Zhang et al. 2019;Tang et al. 2021b).Suppose the same friction parameters are used in this study, ruptures will be challenging to initiate in the designed nucleation zone, because relative strength parameter S (Bizzarri et al. 2011) will be larger.Besides, the steeper dip angle near the surface of the WMF (Fig. 6) tends to increase the normal stress, enhancing the fault clamping effect and making it more difficult to slip.Thus, if our model of the WMF with a steeper dip angle is correct, the WMF will be more stable than the BCF in the long term.

Explanation of rupture patterns
We simulated large earthquake scenarios (models E1-C1).Two interesting scientific issues are worthy of further discussion.First, why do cascading ruptures jump from oblique to strike-slip faults, but the reverse process is challenging?Second, why does the eastward rupture cause greater coseismic displacements than the westward Fig. 13 Evolution of stress and slip for the cases with different nucleation positions.The location of fault point is indicated in Fig. 8 propagation?Although we have mentioned that the rupture directivity can answer the second question, we hope to further explore the deeper reasons from a mechanics perspective.
We answer the first question by calculating the shear stress changes �τ , normal stress change �σ (posi- tive in compression), and the CFS on the two faults under the state when a rupture is approaching the fault intersection (Fig. 12).The CFS is calculated from �CFS = �τ − µ s �σ (King et al. 1994;Harris 1998), measuring whether the fault is brought closer or farther from failure after slip on an adjacent fault (Freed 2005;Parsons et al. 2008;Liu et al. 2018).Ruptures nucleating at the east or west ends of the WMF approach the fault intersection and induce a positive CFS on the LXF.The primary contribution of the increase in CFS comes from �τ rather than �σ , although the normal stress change affects a larger area of the LXF.On the other hand, simulated ruptures nucleating at the north or south ends of the LXF approach the fault intersection and can induce a negative CFS on the WMF.The CFS in most areas near the intersection of the WMF is negative, although there are small-scale positive CFS locally.
The fault activity depends on its prestress state besides the CFS .However, it is more challenging to observa- tionally determine the absolute prestress level or when in its seismic cycle than the principal stress orientation and the stress ratio.The prediction can be more complicated if we consider further stress heterogeneity than those in our modeling.With such a limitation, the LXF has higher initial shear stress than the WMF (Additional file 1: Fig. S4); therefore, the LXF will be closer to failure than the WMF.Consequently, with the evolution of the CFS of the four cases, it is reasonable to expect the cascading ruptures to jump from oblique to strike-slip faults rather than the reverse.
We now answer the second question: Why does the eastward rupture cause greater coseismic slips than the westward propagation?We emphasize that all conditions in the two cases are the same except for the different nucleation positions or rupture propagating directivity.The increase in the coseismic slip is distributed over most of the fault rather than a local area, meaning that the stress evolution at a specific point (the point shown in Fig. 8) on the fault might help reveal the physical mechanism of this phenomenon (Fig. 13).The shear stress variations are similar in the two cases, and the dynamic stress drops are the same.However, the normal stress suddenly increases when the rupture propagates from east to west (Fig. 13a).The seismic waves generated by the rupture superimpose on the edge of the rupture front and couple with the hanging wall moving in the opposite direction (oblique and dextral slip), forming a local compression.This transient increase in normal stress weakened the slip rate to reach a large peak value (Fig. 13c), resulting in a lower coseismic slip than that of reverse propagating directivity.Our results correlate with the study by Tang et al. (2021a).

Conclusions
We developed a new algorithm based on the Wallace-Bott hypothesis for the inversion of fault geometry to infer nonplanar fault geometry in the Wenchuan region.Without sufficient geophysical or geological data to constrain the fault geometry, this is an effective scheme to refine fault geometry with long-term slip-rate data.
Combining deterministic simulations of earthquake dynamics, we assessed the rupture scenarios expected from the present stress field and given the fault geometry of the WMF and LXF.Areas of large ground motion largely depend on where the earthquake nucleates.Several fault nucleation points, friction coefficients, and initial stress states were evaluated, and the general rupture patterns for these earthquake scenarios could fall into three groups.Depending on the initial conditions, the dynamic rupture might start in the LXF, leading to magnitude-7.0earthquakes, or start in the WMF, then cascade through the LXF, leading to magnitude-7.5 earthquakes, or both start and arrest in the WMF, leading to magnitude-6.5 or -7.0 earthquakes.We limited the extent of faults based on available surface fault data, and the maximum magnitudes of future earthquakes can be more significant if faults not in this model interact.
Furthermore, we find the general tendencies of a reverse oblique fault with a conjugate strike-slip fault.Due to the positive CFS near the fault intersection, a rupture starting on the reverse oblique-slip jumps to the strike-slip fault.However, the reverse process is challenging, because the negative CFS dominates near the intersection.Furthermore, the ground motion might be enhanced if the rupture propagation direction is consistent with the movement direction of the hanging wall, because a local decrease of normal stress enhances the coseismic slip.
While we applied our new systematic and physically based method to the case of WMF/LXF, this would apply to any other fault systems in the world where available data are limited to develop models.Therefore, we expect this method to be further tested in other cases and to be able to contribute to the assessment of future earthquake scenarios.(Xu et al. 2009;Liu-Zeng et al. 2009).The blue solid and dotted lines represent the WMF and Lixian faults, both of which did not slip in the main shock but are very active in aftershocks.Accurate fault trace of the WMF is from the field investigation by Xie et al (2011), and the location of the LXF trace is from the distribution of aftershocks (Yin et al. 2018;Li et al. 2019)."Beach ball" represents the lower hemisphere projection of the focal mechanism (strike 225°/dip 39°/rake 120°) of the M w 7.9 ( M s 8.0 ) Wenchuan earthquake obtained from Zhang (2009).(From Li et al. 2019).Figure S6.Inferred regional stress field.a The trade-off between model length and data misfit, e = 1.6 is selected to be the optimal damping parameter.b An example to show the distribution of Closeness between best stress tensor and inverted stress tensors from 1000 times bootstraps in grid number 1.To estimate the stress uncertainty with 95% confidence, we pick 95% of them that have the highest scalar product with the best stress tensor.c Results of crustal stress field inversion from 440 fault plane solutions for the Longmen Shan region with a grid interval of 0.2° both in longitude and latitude (only the stresses in our calculation domain are shown).d The principal stress axes obtained in each block are denoted by red ( σ 1 ), blue ( σ 2 ) and green ( σ 3 ) in the lower hemisphere stereonet.

Fig. 1
Fig.1Map showing the setting of major faults of the center Longmen Shan fault zone.The red bold lines indicate the field-confirmed rupture traces in the Wenchuan earthquake(Xu et al. 2009;Liu-Zeng et al. 2009).The blue solid lines represent the WMF and LXF, both of which did not slip in the mainshock but are very active in aftershocks.The green solid lines represent three active faults parallel to the LXF.Accurate fault trace of the WMF is derived from field investigation ofXie et al (2011), and the location of the LXF trace is deduced from the distribution of pure strike-slip aftershocks(Li et al. 2019).Historical earthquake catalog denoted by filled circles refers toWang et al (2021) and China Earthquake Networks Center (http:// data.earth quake.cn).Green hollow circles denote the historic strong earthquakes (M > 5) before 2008 Wenchuan event.The purple and gray coupled arrows depict the maximum and minimum horizontal compressional stress fields(Luna and Hetland 2013).The third primary stress (not shown) is vertical.The purple dashed regions, respectively, show the positions of Xuelongbao (north) and Pengguan (south) Massif(Shen et al. 2019).The blue box shows the seismic gap near Dayi town.WC-Wenchuan.MX-Maoxian.WMF-Wenchuan Maoxian Fault.LXF-Lixian Fault.BCF-Beichuan Fault.PGF-Pengguan Fault.XYDF-Xiaoyudong Fault.WCEQ-Wenchuan Earthquake.LSEQ-Lushan Earthquake Fig.1Map showing the setting of major faults of the center Longmen Shan fault zone.The red bold lines indicate the field-confirmed rupture traces in the Wenchuan earthquake(Xu et al. 2009;Liu-Zeng et al. 2009).The blue solid lines represent the WMF and LXF, both of which did not slip in the mainshock but are very active in aftershocks.The green solid lines represent three active faults parallel to the LXF.Accurate fault trace of the WMF is derived from field investigation ofXie et al (2011), and the location of the LXF trace is deduced from the distribution of pure strike-slip aftershocks(Li et al. 2019).Historical earthquake catalog denoted by filled circles refers toWang et al (2021) and China Earthquake Networks Center (http:// data.earth quake.cn).Green hollow circles denote the historic strong earthquakes (M > 5) before 2008 Wenchuan event.The purple and gray coupled arrows depict the maximum and minimum horizontal compressional stress fields(Luna and Hetland 2013).The third primary stress (not shown) is vertical.The purple dashed regions, respectively, show the positions of Xuelongbao (north) and Pengguan (south) Massif(Shen et al. 2019).The blue box shows the seismic gap near Dayi town.WC-Wenchuan.MX-Maoxian.WMF-Wenchuan Maoxian Fault.LXF-Lixian Fault.BCF-Beichuan Fault.PGF-Pengguan Fault.XYDF-Xiaoyudong Fault.WCEQ-Wenchuan Earthquake.LSEQ-Lushan Earthquake

Fig. 4
Fig. 4 Schematic to describe how the fault geometry is related to traction at measured points.We get the strike and dip angles for each collocation point using cubic spline interpolation (McKinley and Levine 1998), and the traction on each collocation point can be calculated through resolving the heterogeneous tectonic stress tensor along the fault strike and dip.Then, we get the averaged traction from all collocation points in the spherical computational region with radius R

Fig. 5
Fig. 5 Map showing the measured points of long-term slip rate along the WMF.Each gray box contains the locations and long-term slip rates measured from different researchers.Dextral means dextral strike-slip rates and Up (NW) means vertical slip rates with northwestern block being hanging wall.The red and blue lines are the same with Fig. 1.The red stars denote the measured points of long-term slip rate.Gray circles represent aftershocks of the 2008 earthquake(Yin et al. 2018)

Fig. 7
Fig. 7 Uncertainty of WMF geometry estimation.a Observed long-term slip rates and scaled traction.The blue and orange lines present weighted average shear stress along the dip and strike, respectively.The length of the bar indicates the standard deviation of slip rate.The vertical thin pink lines indicate the position of B-spline curves.b Distribution of inverted parameters τ, the red line is the average of 150 inversions (gray lines), and blue triangles indicate the position of interpolation curves.The average of inverted fault depth h 0 = 16.11 is not shown here.c Objective functions versus iterations from 100 inversions.d Standard variance distribution of fault dips derived from 150 inversions

Fig. 8
Fig. 8 Total slips and earthquake magnitudes from different nucleation points noted by gray arrows.The green arrow denotes the point in Fig. 13

Fig. 9
Fig. 9 Simulations (model E3 ~ C3 in Table 1) of the case 2 for depth-dependent vertical stress, the lower right sub-figure shows two different depth-dependent stress configurations with saturation depth 5 km and 10 km

Fig. 10
Fig. 10 Slip patterns from all models using simple schematic diagrams a and histogram of the magnitudes ( M w ) obtained in the 20 simulations b.Yellow highlighting indicates the path of propagating rupture, and solid red circles denote the nucleation sites
FigureS3.Boundary elements for each fault segment and free surface.We employed the boundary elements of approximately 26,824, and control triangular elements to be 500 m on the fault and free surface in the vicinity of the fault, and increase this size with a maximum rate of 7.5 times toward the outer boundary of the surface until it reaches the maximum size of around 2000 m. Figure S4.Initial on-fault shear stress and normal stress resolved by the tectonic stress, as well as Initial on-fault shear stress drop τ drop = τ 0 − µ d σ , μ d = 0.12. Figure S5. a The spatial distribution of 391 focal mechanism solutions (beachball, lower hemisphere stereographic projection) determined in this study from earthquakes between 2009 and 2016.The color and size of the beachballs indicate the type and seismic moment, respectively.Red, normal (NF) or normal oblique faulting (NS); green, strikeslip faulting (SS); blue, thrust (TF) or thrust oblique faulting (TS); black, mixed type.b Histogram illustrating the faulting type distribution of the earthquake focal mechanism data determined in this study.

Figure S7 .
Along strike distribution of scaled three-component principal stress.The scaled three-component principal stress are marked red ( σ 1 ), blue ( σ 2 ) and green ( σ 3 ) colors, which are derived using bootstrap resampling of the selected fault planes.The black circles denote the stress component closest to the vertical, and the dotted lines are the average of each stress component (a).The orange solid line and red dash line represent the trend and average of principal stress orientation (b).

Figure S8 .
Inversions of fault geometry with different horizontal intervals of interpolated curves.From top to bottom: Observed long-term slip rates and scaled traction; Distribution of inverted parameters η ; Objective func- tions versus iterations; Final fault surface derived from inverted parameters η .Other details for each plot refer to Fig. 4.

Figure S9 .
Same as Additional file 1: Fig.S7, but with different integral radius.
Figure S11.a Calculated Coulomb failure stress change and the movement of Lixian Fault, assuming µ s = 0.52 (Modified from Tang et al. 2021a, b).The spatial distribution of aftershocks is from Yin et al. (2018).b The Coulomb stress change of the WMF, assuming that the LXF ruptures at southern end.Figure S12.Three components and magnitude of surface slip, with eastern and western nucleation points.The white pentagrams denote the epicenter.The white arrow indicates the NE direction.See Fig. 5 for the definition of the x-y-z components.
Figure S13.Magnitude of peak ground velocities from different nucleation points.The white pentagrams denote the epicenter.

Figure A1 .
Schematic to describe the calculation of traction from interpolated fault surface.Green circles are the control points at equal intervals on the dip curves, and orange circles denote the equally spaced points on the interpolation curves.The average interval of collocation points is around 3 km.The inverted triangles represent the locations of the measured points.

Additional file 1: Text S1.
Stress inversion.Text S2.Test on parameters for fault geometry inversion.Text S3.The effect from Wenchuan mainshock .Appendix A. Conjugate gradient method for fault geometry inversion.Text S4.Peak ground displacement and velocity.Figure S1.Map showing the setting of major faults of the center Longmen Shan fault zone.The red bold black lines indicate the field-confirmed rupture traces in Wenchuan earthquake