Constraining Shear Strength of Fault Damage Zone Using Geodetic Data and Numerical Simulation

Shear strength of damage zone, representing the stress threshold for rupture initiation, is a critical parameter in faulting mechanics. Despite its significance, the damage‐zone's shear strength has not been estimated in natural earthquake ruptures. Here we employed coseismic deformation and strain, kinematic slip model, and finite element modeling to determine the elastic properties and peak shear stress of coseismic damage zones along the 2021 Mw 7.4 Maduo earthquake. Through the analysis of the lowest shear stress resulting in surface ruptures and the highest stress without surface rupture, we constrained the strength within a range of 7–17 MPa. Our result is consistent with strength (5–16 MPa) of sandstone samples from laboratory tests, demonstrating the validity of this estimation. Although factors such as fault maturity and confining pressure influence strength variation, the strength can directly reflect the stress threshold required for macroscopic surface rupture formation in fault damage zones dominated by sandstone.


Introduction
For a damage zone embedded in a seismogenic fault, the shear strength is one of its most significant mechanical properties.The shear strength represents the ability that a damage zone can withstand maximum shear stress without generating macroscopic ruptures (Bletery et al., 2016;Faulkner et al., 2010;Fossen, 2016;Fossen & Cavalcante, 2017;Marone et al., 1995;Tenthorey et al., 2003).Compared to the surrounding host rocks, damage zones often exhibit a dense network of fractures, which can be macroscopically approximated as an elastic zone with a reduced elastic modulus (e.g., shear and tensile strength, and stiffness) and seismic velocities (Allam & Ben-Zion, 2012;Biegel & Sammis, 2004;Chester & Logan, 1986;Choi et al., 2016).Damage zones have been found to cause a notable 20%-60% decrease in seismic velocities (Lewis & Ben-Zion, 2010;Lu & Ben-Zion, 2022;Perrin et al., 2016), and to accommodate 28%-88% of distributed, off-fault deformation (Antoine et al., 2021;Dolan & Haravitch, 2014;C. Li et al., 2022;Zinke et al., 2014).Such mechanical response of damage zones during the earthquake cycle may in principle be related to mechanical properties anomalies within the damage zone (Ishii, 2015;Okubo et al., 2019;Shearer, 2019;Thakur et al., 2020).Therefore, understanding shear strength of a damage zone holds significant importance for variable applications, such as assessing rupture mechanics and hazards associated with shallow faults.
The prevalent method to constrain the mechanical properties of rock materials is to conduct laboratory experiments (Hart & Wang, 1995;Ishii, 2015Ishii, , 2016;;Wong & Baud, 2012;Xu, Fukuyama, et al., 2023;Zhou et al., 2018).However, even after a range of tests on host rocks with different lithologies, confining pressures, and fluid saturations, the shear strength derived from small-scale, intact host rock samples in the laboratory often fails to accurately represent that of larger rock masses (Jónsson, 2012).Other methods for estimating the mechanical properties of damage zones include numerical models (e.g., Bletery et al., 2017;Fang & Dunham, 2013;Huang, 2018;Thakur et al., 2020) and seismic dense arrays (e.g., Zhang et al., 2022;Zhou et al., 2022).Importantly, the shear strength of a damage zone can be significantly influenced by many factors, such as fault maturity (Ikari et al., 2011;Manighetti et al., 2007;Thakur & Huang, 2021), thermal dependence of fault strength (Badt et al., 2020;Di Toro et al., 2011), fault surface topography (Xu, Fukuyama, et al., 2023), and dynamic weakening law (Fang & Dunham, 2013;Xu et al., 2018).Existing methods fail to replicate structures of natural damage zones in a laboratory setting and do not account for dynamic rupture effects during natural earthquakes.Therefore, the assessment of "how much shear stress a fault damage zone can withstand before generating macroscopic rupture" is hindered by the lack of direct estimates from natural earthquake ruptures.
Notably, recent advancements in satellite observations for coseismic deformation and strain have significantly improved our understanding of mechanical properties within fault damage zones, such as fault friction, inelastic failure, rupturing strain limit, and elastic modulus (e.g., Barnhart et al., 2020;C. Li et al., 2023;Milliner et al., 2022;Xu et al., 2020;Xu, Liu, & Lavier, 2023).These efforts enable us to characterize the coseismic strain of the damage zone in great detail, a crucial prerequisite for resolving damage-zone's shear strength in natural earthquake ruptures.Another prerequisite demands the presence of both surface-rupture segments and segments without surface rupture (referred as non-rupture) during a single earthquake.The overlapping range of coseismic shear stress observed along the two types of segments affords us a unique opportunity to directly determine the shear strength of the damage zone.Interestingly, this specific requirement appears to have been met in the 2021 Maduo earthquake (eastern Tibetan Plateau).The Maduo seismogenic fault produced surface rupture segments interspersed with non-rupture segments of a similar spatial scale (Liu-Zeng et al., 2022;Pan et al., 2022;Yuan et al., 2022).In our study, we employed the coseismic deformation and strain, kinematic slip model, and finite element modeling to determine the structural parameters and elastic properties of coseismic damage zones.By combining the maximum shear strain along the Maduo rupture (C.Li et al., 2023), we quantified the peak shear stress imposed by the Maduo rupture, and directly constrained the shear strength of the damage zone.

Methods of Constraining Shear Strength
Rock deformation commonly undergoes distinct stages defined by two thresholds: yield point and rupture point (Fossen, 2016;Jónsson, 2012; Figure 1).The yield point (i.e., yield strength) marks the transition from elastic to plastic deformation.The rupture point indicates the stress limit the rock can withstand before forming macroscopically fracturing, with the shear stress limit aligning with the shear strength.For the Maduo earthquake, the occurrence of surface rupture depends on whether the peak coseismic shear stress accommodated within the damage zone surpasses its shear strength (Figure 1a).Along the rupture segment R1-6, the coseismic shear stress (σ rupture ) has exceeded both the yield and rupture point of the damage zone because of the presence of surface ruptures (Figure 1b).In contrast, along the non-rupture segment N1-6, the coseismic shear stress (σ non-rupture ) has not exceeded the rupture threshold of the damage zone because there were no surface ruptures, although it may have exceeded the yield threshold (Figure 1c).We can approximately determine the peak coseismic shear stress using the maximum cosiesmic shear strain (ε non-rupture and ε rupture ) and elastic shear modulus (G) based on quasielastic stress-strain relationship (Figures 1b and 1c).
However, strain hardening results in the occurrence of some inelastic deformation that has a non-linear stressstrain relationship between the yield and rupture points.The previously-mentioned approximation may result in an overestimation of the peak shear stress.However, this overestimation only has limited impact on our stress results.The rationale behind is we further consider the following two key factors.Firstly, near-surface bedrocks primarily undergo brittle deformation mechanism, typically experiencing brittle ruptures once their yielding threshold is exceeded, with little to negligible plastic strain (Cartwright-Taylor et al., 2022;Meyers & Chawla, 2008;Nouri et al., 2022;Wang et al., 2020;Wong & Baud, 2012;Zhou et al., 2018).This implies that the quasi-elastic stress-strain relationship is applicable before rupture occurs.Secondly, the coseismic shear strain is continuous along the entire fault, however they are localized and highly distributed along the rupture and nonrupture segments, respectively (Antoine et al., 2023;C. Li et al., 2023;Figure S1 in Supporting Information S1).As a result, the distributed shear strain of non-rupture segments is insufficient to trigger brittle surface

Writing -review & editing: Xinjian Shan
Geophysical Research Letters 10.1029/2024GL108169 rupture (Liu-Zeng et al., 2022;Pan et al., 2022;Yuan et al., 2022).While some plastic components might be present, coseismic displacements in the non-rupture segments are more likely to be dominated by the quasi-elastic regime.Consequently, the shear stress estimated from the quasi-elastic equation can effectively represent the peak value imposed by the Maduo rupture.
After determining the peak coseismic shear stress, the lowest σ rupture in each rupture segment represents the minimum coseismic shear stress required to initiate surface rupture (Figures 1b and 1d).Conversely, the highest σ rupture in each non-rupture segment represents the maximum shear stress that damage zone can withstand without developing surface rupture (Figures 1c and 1e).The overlapping range of the lowest σ rupture and the highest σ non- rupture represents a critical state of shear stress.In the cases where damage zone experiences coseismic shear stress above this range, surface rupture occurs; whereas surface rupture does not occur.Therefore, the overlapping range of the lowest σ rupture and the highest σ non-rupture directly can indicate the shear strength of coseismic damage zone.

Estimating Elastic Properties of the Damage Zone
The estimation of peak shear stress requires quantifying both the maximum shear strain and shear modulus of coseismic damage zones.The maximum shear strain along the Maduo rupture has been measured from the shear strain fields derived from the SPOT satellite optical images (C.Li et al., 2023).Therefore, we only need to determine the shear modulus of coseismic damage zones.To this end, we utilized the geodetic finite element modeling approach from a prior study of Xu, Liu, & Lavier, 2023.Do note that with an assumption of an isotropic linear elastic medium, the system of governing equations is given by linear elasticity.We selected the coseismic displacements in the non-rupture segments N1-6 as our modeling target because these displacements are more likely to be dominated by the elastic regime, therefore more suitable for elastic modeling.
To determine the structural parameters of coseismic damage zones across the segments N1-6, we used five acrossfault profiles AA′-EE to extract the fault-parallel displacement from the optical-based displacement field (Figure 2a).The location and width of each coseismic damage zones in the segments N1-6 can be measured from five displacement profiles (Figure 2b).Additionally, to estimate dip angle and depth of each damage zone, we constructed the fault slip model of the Maduo earthquake using the three-dimensional coseismic deformation resolved from both Sentinel-1 and ALOS-2 images (Liu et al., 2022; Text S1 in Supporting Information S1; Figure 2c and Figures S2-S4 in Supporting Information S1).We determined the structural parameters of each damage zones, including the locations and widths along the profiles AA′-EE′, as well as corresponding depths and optimal dip angles at the fault segment 1, 4, 8, 11, and 13 in our kinematic slip model (Figure 2d).
Based on the structural parameters of five damage zones, we further discretized the meshes using an open-source mesh generator Gmsh (Geuzaine & Remacle, 2009).The target damage zones were embedded within a virtual three-dimensional host rock with size of 8 km in the x-, y-, and z-directions, respectively.The shear modulus of the damage zone was then estimated through determining shear modulus reduction ratio (G d /G h ) between intact host rock with a shear modulus G h and damage zone (Xu, Liu, & Lavier, 2023).We denote an isotropic linear elastic host rock medium as V 1 -V 3 .The first scenario is that a shear stress was subjected to an intact rock V 4 embedded within the host rock V 1 -V 3 (Figure S5a in Supporting Information S1).The second scenario is that the same shear stress was imposed to a damage zone V 4 embedded within the same host rock V 1 -V 3 (Figure S5b in Supporting Information S1).Such differential amount of the shear displacement in the two scenarios should be contributed to reduction in shear modulus between intact host rock and damage zone (Figure S5c in Supporting Information S1).We implemented the FEM to simulate shear displacement in FEniCS, an open-source computing platform for solving partial differential equations (Alnaes et al., 2015).The boundary conditions of each model ensured that the strike-slip and dip-slip components matched the estimates from our slip kinematic model (Figure S5 in Supporting Information S1).The bulk and shear modulus for the isotropic materials can be readily wellsolved, allowing us to derive the Poisson's ratio (v) from the isotropic stress-strain relationship.To ensure the accuracy of the computational results and reduce computation time, we used an iterative solver based on the Krylov method, with a convergence tolerance of 10 14 .
Given that the Maduo earthquake rupture passed through a low-relief mountainous area where Permian-Triassic sandstone is mostly exposed at the surface (Cheng et al., 2006;C. Li et al., 2023; Figure 1), we adopted a Young's modulus of 15 GPa for an intact sandstone based on laboratory experiments, and its shear modulus was further estimated to be 6 GPa based on Poisson's ratio of 0.25 (Fossen, 2016;Nouri et al., 2022).Next, we searched for the optimal G d /G h and v, which can result in the best match between the shear displacement simulated by the FEM and the observed surface displacement for the across-fault profiles AA′-EE′ using a MCMC Bayesian sampler (Foreman-Mackey et al., 2013).The goodness of fit between our observed surface displacement and FEM displacement was evaluated using the normalized chi-square misfit.We assumed a uniform prior probability distribution for G d /G h and v. Our Monte Carlo Chain produced 1,000 samples of the postprior distribution.The points with highest probability density were used to determine the optimal G d /G h and v.

Results
The marginal posteriors of shear modulus reduction ratio G d /G h and Poisson's ratio v are shown in Figure 3a and Figure S7a in Supporting Information S1.After exploring a range of model parameters, we found that the five preferred shear modulus reduction ratio is 0.10-0.27,and the corresponding Poisson's ratio is 0.25-0.30(Figure 3b and Figure S7b in Supporting Information S1).This suggests that the strengths of the coseismic damage zone have been decreased to 10%-27% of the intact rock strength.Since the shear modulus of the intact sandstone is assumed to be 15 GPa, we could constrain the shear modulus G d of the damage zone with a range of 0.6-1.6GPa (Figure 4a).These shear displacements simulated by the FEM using the preferred model parameters excellently fitted the observed coseismic displacement (Figure 3c and Figure S7c in Supporting Information S1).Subsequently, we collected the maximum shear strain (ε rupture and ε non-rupture ) along the Maduo rupture (C.Li et al., 2023;Figure 4a).Together with our inferred shear modulus of each damage zone, we further estimated the peak coseismic shear stress along the rupture segments (σ rupture = G d ε rupture ) and the non-rupture segments (σ non-rupture = G d ε non-rupture ) imposed by the Maduo rupture.
For each rupture segment, the lowest σ rupture imposed by the Maduo coseismic rupture was determined to be the minimum value among all σ rupture , with a range of 7-15 MPa (Figure 4b).The lowest σ rupture represents the minimum shear stress required to trigger surface rupture within the damage zone.In contrast, all σ non-rupture along non-rupture segments yield six highest values ranging from 8 to 17 MPa (Figure 4b).The highest σ non-rupture signifies the capacity of the damage zone to withstand the utmost shear stress without initiating surface rupture.The overlapping range of the lowest σ rupture and highest σ non-rupture in the contiguous pairs directly indicates the shear strength of the damage zone.Therefore, the shear strength of the damage zone varies; it ranges from approximately 14 to 17 MPa near the western (R1-N1-R2), and 12-13 MPa near eastern rupture ends (R6-N6).In the segments R2-N3, the shear strength varies between 7 and 11 MPa.Similarly, from segments N3 to R5, four bounds indicate shear strength values ranging from 8 to 17 MPa.In the segments R5-N5, the shear strength can be constrained to range from 7 to 13 MPa.In general, based on the joint constraints from the five contiguous pairs of damage zones, the shear strength required for the formation of coseismic surface rupture should fall within a range of 7-17 MPa.Note that our given shear strength is relative, with the range of 7-17 MPa being derived based on the premise of a Young's modulus of 15 GPa for intact sandstone.

Discussions and Conclusions
Our study provides a direct estimate of the shear strength for fault damage zones based on a natural earthquake rupture.Notably, the geophysical significance of our strength represents the maximum shear stress that the damage zone can withstand before generating macroscopic rupture at surface.However, it is essential to note that that our strength estimated in this study only represents static strength at surface or shallowest surface, as surface bedrocks are likely to be significantly lower rigid than rocks at several kilometers deep in the crust (Fossen, 2016;Jónsson, 2012;Rubin & Pollard, 1987).Given the proportionality of shear strength to crust depth, our finding can be further applied to investigate the faulting strength at seismogenic depth (e.g., 5-15 km).Additionally, the Maduo rupture predominantly traversed the Permian-Triassic sandstone, exhibiting no significant along-strike lithological variations (Cheng et al., 2006;C. Li et al., 2023; Figure 1).The damage-zone's shear modulus (G d ) and strength are estimated using the 15 GPa Young's modulus of host sandstone (G h ).Given similar Young's modulus among rocks of the same lithology under normal temperature and pressure, employing a single Young's modulus to estimate shear strength along the Maduo rupture is justifiable.However, we emphasize that our inferred strength values should still be relevant for damage zones that have developed in sandstone-dominated areas.For the damage zone embedded in other rock lithology, it is necessary to use the corresponding Young's modulus to re-estimate shear strength.We found the shear modulus of the damage zones along the Maduo rupture shows an overall increasing trend from the western to the eastern rupture end (Figure 4a).Previous studies have illustrated that the mechanical properties of fault damage zones systematically depend on fault maturity; young developing faults may exhibit relatively high strength, whereas mature well-slipped faults can be weaker, possibly due to activation of various weakening mechanisms with an increasing cumulative slip (Collettini et al., 2009;Fialko, 2021;Fossen, 2016;Ikari et al., 2011;Manighetti et al., 2007;Thakur & Huang, 2021).It is noted that the Kunlunshankou-Jiangcuo Fault (i.e., the seismogenic fault of the Maduo earthquake) connects with the Eastern Kunlun Fault at its western end and terminates at its eastern end (coincident with the Maduo eastern rupture end), where the slip accumulation is negligible (C.Li et al., 2023;Pan et al., 2022).Consequently, this fault likely has developed through progressive eastward propagation, and therefore, its slip accumulation and structural maturity should decrease toward the east (e.g., Manighetti et al., 2001;Perrin et al., 2016).Given that the rock within the fault damage zone tends to weaken, as fault slip accumulates (e.g., Ben-Zion & Sammis, 2003;Chester & Logan, 1986), the increase in shear modulus toward the east can be explained by the decrease in fault maturity in that direction.The fault zone closer to the eastern end is expected to have lower slip accumulation and maturity, resulting in a less developed gouge zone and, consequently, higher shear strength.
Apart from the fault maturity, the along-strike variation of the shear strength may be due to several potential factors, including rock lithology, coseismic dynamic weakening, and loading conditions like confining pressure, temperature, and the presence of fluids (Figure 4b; Badt et al., 2020;Di Toro et al., 2011;Fang & Dunham, 2013;Jamtveit et al., 2019;Wang et al., 2020;Xu, Fukuyama, et al., 2023;Xu et al., 2018).Along the Maduo earthquake fault, the predominant bedrock consists of Permian-Triassic sandstone (Figure 1).Despite the presence of Quaternary deposits covering on the bedrock, our strength values should still be relevant for damage zones that have developed in sandstone-dominated areas.For example, the strength (8-15 MPa) from the bedrock terrain between the segments N4-R5 is seemly lower than that (12-13 MPa) from the Quaternary alluvial terrain at the segments R6-N6.This contradicts the expectation that the terrain with softer Quaternary deposits has a lower strength (Figure 4; Cheng et al., 2006).This variation of shear strength cannot be attributed solely to the difference in rock lithology.
The confining pressure has a significant sensitivity to the rock mechanical properties.With the increase of confining pressure, the rock deformation undergoes transformation process from brittle failure to ductile deformation, and the shear strength also increases (Wang et al., 2020).At a rough estimate, the Quaternary deposits with a thickness of <200 m would impose crustal confining pressure of <5 MPa (lithostatic stress σ = ρgh = 2,584 kg/m 3 × 9.8 m/s 2 × 200 m) to the coseismic damage zone (Cheng et al., 2006;Fialko, 2021;Huang et al., 2020;Figure 1a).For example, (a) the parts R2-N2 covered by the Quaternary deposits was inferred a strength range of 8-11 MPa, and nearby bedrock-terrain segments R3-N3 has a relative lower strength of 7-8 MPa; (b) the strength of the rupture segments (i.e., R1, R2, and R6); covered by the Quaternary deposits is significantly larger than that at the bedrock-terrain rupture segments (Figure 4b).This suggests that the confining pressure imposed by the Quaternary deposits positively influences the shear strength.The relatively low confining pressure of 5 MPa does typically not change the rock brittleness such that our above assumption is still applicable (Wang et al., 2020); however, as confining pressure increases, the shear strength of the damage zone is also likely to rise.
Except for the factors mentioned, temperature can also have an impact on fault strength (Fossen, 2016).In our investigation, the near-surface areas typically remain at ordinary temperature conditions, and the temperature changes are not significant enough to cause substantial variations in fault shear strength.However, it is necessary to take into account coseismic dynamic weakening that can occur due to high slip rates (Chang et al., 2012;Xu et al., 2018).For example, at high seismic slip rates, there is a remarkable decrease in the friction coefficient of crustal silicate rocks due to intense "flash" heating of microscopic asperity contacts, leading to a degradation of their shear strengths (Goldsby & Tullis, 2011).During the Maduo rupture, the estimated coseismic slip rate was approximately (at most) 0.4 m/s near the E30 km, decreasing to less than 0.1 m/s at both the western and eastern rupture ends.Paradoxically, the strengths at these ends are higher, which is in clear contradiction to the expected degradation of fault shear strength caused by a fast slip rate.The variation in strength does not appear to be attributable to the seismic slip rate either, as the inverted coseismic slip rates are remarkably consistent (0.2-0.3 m/s) across all measurement sites (Fang et al., 2022).Furthermore, our derived shear strengths are comparable to laboratory measurements on the intact siliceous mudstones and Cenozoic-Mesozoic sandstones.When considering shear modulus reduction of damage zone compared to intact rocks, the shear strength of 7-17 MPa is Geophysical Research Letters 10.1029/2024GL108169 consistent with those of intact siliceous mudstone samples (6-19 MPa) and sandstone samples (5-16 MPa) tested on various conditions (Ishii, 2015;Nouri et al., 2022;Figure S8 in Supporting Information S1).Most importantly, compared to laboratory measurements, our estimations provide a more comprehensive reflection of natural damage zones' ability to withstand coseismic shear stress, as they incorporate interactions of various weakening mechanisms and tectonic conditions during natural earthquake ruptures.
In conclusion, we employed coseismic deformation and strain, kinematic slip model, and finite element modeling to determine the elastic properties and peak shear stress of coseismic damage zones.Our results revealed that for the damage zone embedded in an intact sandstone with a Young's modulus of 15 GPa, its shear strength can be estimated to range from 7 to 17 MPa.Our findings also suggest that shear strength of damage zone is positively correlated with confining pressure and negatively with fault maturity, respectively.Overall, our estimations reflect the capacity of natural damage zones in sandstone-dominated areas to withstand coseismic shear stress without initiating macroscopic rupture.

Figure 1 .
Figure 1.(a) Coseismic surface rupture (R1-6) and non-rupture segments (N1-6) generated by the Maduo earthquake, with geological units.Modified from C. Li et al. (2023).(b, c) Conceptualized shear stress-strain curves of coseismic surface rupture and non-rupture segments, respectively.(d, e) Schematics of coseismic damage zone for surface rupture and non-rupture segments, respectively.Notably, surface rupture can occur when peak coseismic shear stress exceeds strength of a damage zone in panel (d); otherwise, it will not occur in panel (e).

Figure 2 .
Figure 2. (a) Coseismic fault-parallel displacement and shear strain fields collected from C. Li et al. (2023).(b) Fault-parallel surface displacement along the non-rupture segments N1-6 extracted from the across-fault profiles AA′-EE′ (locations shown in panel a).(c) Coseismic slip distribution of the Maduo fault.(d) FEM setup for analysis of the surface slips across the profile AA′-EE′ buried in a homogeneous medium.The sketch of the model meshes was generated by Gmsh, with solution computed by FEniCS and visualized in PyVista (Sullivan & Kaszynski, 2019).

Figure 3 .
Figure 3. (a) Marginal posteriors of shear modulus ratio (G d /G h ) and Poisson's ratio (v) for the damage zones along the profile AA′-CC′.The red dots represent the best matching v and G d /G h (±1σ).(b) Displacement distribution in each FEM using the best matching parameters.Remarkably, the inferred shear modulus totally shows an increasing tendency toward the east.(c) Observed and FEM displacement profiles.The black and red lines respectively denote the amplitude of shear displacement for the input model and the best matching samples.Notably, the optimal simulated displacement closely aligns with the observed data.

Figure 4 .
Figure 4. (a) The maximum shear strain along the Maduo earthquake rupture, and the shear modulus of the damage zones.(b) The estimated peak coseismic shear stress imposed by the Maduo rupture.The shear strengths of the damage zones were inferred from the lowest shear stress at the surface rupture segments and the highest at the non-rupture segments.