# On the importance of 3D stress state in 2D earthquake rupture simulations with off-fault deformation

Louise Jeandet, Marion Thomas, Harsha Bhat

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

## To cite this version:

Louise Jeandet, Marion Thomas, Harsha Bhat. On the importance of 3D stress state in 2D earthquake rupture simulations with off-fault deformation. EGU General Assembly 2023, European Geosciences Union, Apr 2023, Vienna, Austria. 10.5194/egusphere-egu23-6708. hal-04189796

HAL Id: hal-04189796

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

Submitted on 22 Jul 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# On the importance of 3D stress state in 2D earthquake rupture simulations with off-fault deformation

Manuscript accepted to be published in Geophys. J. Int.

Louise Jeandet Ribes¹, Marion Y. Thomas¹ and Harsha S. Bhat²

1. Institut des Sciences de la Terre Paris, Sorbonne Université, CNRS-UMR 7193, Paris, France.

2. Laboratoire de Géologie, École Normale Supérieure, CNRS-UMR 8538, PSL Research University, Paris, France.

During the last decades, many numerical models have been developed to investigate the conditions for seismic and aseismic slip. Those models explore the behavior of frictional faults, embedded in either elastic or inelastic media, and submitted to a far field loading (seismic cycle models), or initial stresses (single dynamic rupture models). Those initial conditions impact both on-fault and off-fault dynamics. Because of the sparsity of direct measurements of fault stresses, modelers have to make assumptions about these initial conditions. To this day, Anderson's theory is the only framework that can be used to link fault generation and reactivation to the three-dimensional stress field. In this work we look at the role of the three dimensional stress field in modelling a 2D strike-slip fault under plane-strain conditions. We show that setting up an incorrect initial stress field, based on Anderson's theory, can lead to underestimation of the damage zone width by up to a factor of six, for the studied cases. Moreover, because of the interactions between fault slip and off-fault deformation, initial stress field influences the rupture propagation. Our study emphasizes the need to set up the correct initial 3D stress field, even in 2D numerical simulations.

**KEYWORDS:** Earthquake dynamics, Dynamics and mechanics of faulting, Elasticity and anelasticity, Numerical modelling, Friction

## 1 Introduction: Modelling fault slip

Catastrophic, seismic events of large magnitude ($M_w > 7$) remain relatively rare, with a recurrence time of several decades up to a millennium (Cubas et al., 2015). As a consequence, observations are sparse and numerical models are powerful tools to explore the conditions for seismic and aseismic fault slip. Numerical models that account for both seismic slip and long-term slow slip are challenging because of the wide range of temporal and spatial scales involved (Lapusta et al., 2000). Hence, compromises are made to reduce the computational cost depending on the goal of the model. These compromises can be safely categorized as models that focus on a single dynamic rupture event and ones that model multiple seismic cycles.

Single dynamic rupture models were some of the first analytical and numerical models developed for earthquakes. They reproduce an event from the moment it turns dynamic to its arrest and have provided important insights into earthquake mechanics (Kostrov, 1964; Andrews, 1976; Madariaga,
1976; Das & Kostrov, 1986, 1987, among many others that followed).

For seismic cycle models, the focus is on integrating all stages of fault slip (inter-, co- and post-seismic) over thousands of years (Erickson et al., 2020). This is indeed critical as pre-stress inherited from aseismic slip and prior seismic events likely play a determinant role on where earthquakes will nucleate and how far their rupture will propagate (Thomas et al., 2014). The most popular strategy is to model a single planar fault governed by rate-and-state (R&S) friction law, embedded in a purely elastic medium. (e.g., Ben-Zion & Rice, 1997; Lapusta et al., 2000; Richards-Dinger & Shearer, 2000; Kato, 2004; Barbot et al., 2012; Im et al., 2020; Liu et al., 2020). In these models, frictional heterogeneities are invariably advocated to reproduce the full slip spectrum (creep, slow slip events, earthquakes, etc...). More recent models also account for complex fault geometry (Romanet et al., 2018; Ozawa & Ando, 2021; Uphoff et al., 2022, among others). However, due to timescales that vary over several decades of orders of magnitude, the most popular compromise that is made is to ignore the wave mediated stress transfer (inertial dynamics) and only account for linear elastic static stress transfer along with the instantaneous, local traction contribution. Such class of models are called quasi-dynamic models.

Natural fault zones, however, i.e. the fault plane and its surrounding medium, are intricate structures in constant evolution in response to tectonic loading. As an example, during an earthquake, part of the stored elastic strain energy is dissipated in off-fault deformation, or damage, which in turn radiate and affect the slip dynamics of the main fault. If we provide an idealized description, fault zones comprise of a non-planar fault core, where most of the displacement has occurred, surrounded by a damage zone that has a spatial scale of the order of meters to kilometers (e.g., Chester et al., 1993; Biegel & Sammis, 2004; Faulkner et al., 2011). Frequently, the fault core corresponds to an extremely narrow band and the damaged wall rock includes layers of gouge and breccia bordered by fractured rocks. The last two layers are included in the damage zone because they lacked extensive shearing. These structures are of key importance in the mechanics of faulting. For example, fault roughness has an effect on the fault resistance to slip, i.e., the fault strength (e.g., Dunham et al., 2011a; Tal & Faulkner, 2022). Laboratory experiments of earthquakes in a damaged medium show that there is an intimate interaction between the rupture and off-fault damage zone (Sammis et al., 2009; Bhat et al., 2010; Biegel et al., 2010). The density of this damage has a direct impact on the elastic moduli of the bulk (Walsh, 1965b,a; Faulkner et al., 2006), therefore, on the quantity of strain energy which is stored and further released by fault slip. In fact, systematic micro- and macrostructural field studies have been performed on damage zones (e.g., Shipton & Cowie, 2001; Manighetti et al., 2001, 2004; Dor et al., 2006; Mitchell & Faulkner, 2009; Savage & Brodsky, 2011; Johnson et al., 2021; Rodriguez Padilla et al., 2022), as a key component to understand the energy balance of earthquakes (e.g., Rice, 2002; Kanamori, 2006; Okubo et al., 2019). Off fault damage is observed from the mesososcopic scale to the microscopic scale, with a microfracture density that decreases exponentially away from the fault core (Mitchell & Faulkner, 2009). The width of the damage zone is also believed to be decreasing with depth, forming a “flower-like structure” (Ben-Zion et al., 2003; Cochran et al., 2009). However Okubo et al. (2019) have demonstrated numerically that, despite its reduction in spatial-extent with depth, energetically speaking, the contribution of the off-fault damage increases with depth.

Thus, as oppose to seismic cycle codes, a second set of models have been developed to catch the full slip dynamics, the wave propagation and the interactions with the off-fault medium during an earthquake. With these models, researchers have explored the effect of complex fault geometry and/or
the effect of off-fault damage on seismic rupture. Some models treat the bulk as a solid linear-elastic material and prescribe a low-velocity zone around the fault to account for damage (e.g., Cappa et al., 2014; Huang et al., 2014). Another set of models has explored the effect of spontaneous dynamic generation of off-fault deformation using a Mohr-Coulomb (e.g. Andrews, 2005; Ben-Zion & Shi, 2005; Hok et al., 2010; Gabriel et al., 2013) or Drucker-Prager (e.g. Templeton & Rice, 2008; Ma, 2008; Dunham et al., 2011a; Johri et al., 2014) based plastic constitutive laws. Another class of models have treated off-fault damage as tensile cracks, using a stress- (Yamashita, 2000) or fracture-energy-based (Dalguer et al., 2003) criterion. Okubo et al. (2019) used a Finite Discrete Element Method (FDEM) that allows spontaneous nucleation and propagation of off-fault fracture network in a medium. These studies have provided a good insight on the effect of a fault zone structure on dynamic ruptures but the models do not account for the observed coseismic change of elastic properties in the bulk (Hiramatsu et al., 2005; Brenguier et al., 2008; Froment et al., 2014) which also impacts the rupture dynamics. This can be achieved by using a homogenized damage mechanics based constitutive law (e.g., Lyakhovsky et al., 1997; Bhat et al., 2012a; Xu et al., 2014; Thomas et al., 2017b; Thomas & Bhat, 2018). A vast majority of the investigations cited above are done under two-dimensional plane-strain conditions.

Finally, some models have been developed to look at the zeroth order effect of the fault core (as opposed to a fault interface) or the damage zone on the seismic cycle. To overcome the computational cost, they have to compromise on both the fault slip dynamics and the dynamics of bulk evolution (off-fault crack growth). As an example, van den Ende et al. (2018) have represented the fault core by a shear band but the bulk remains elastic. Kaneko et al. (2011); Idini & Ampuero (2020), and Abdelmeguid et al. (2019) among others, looked at the effect of a prescribed low velocity fault zones (LVFZ), but, by construction, the bulk is still a passive elastic body. They used quasi-static or quasi-dynamic approximations to solve the problem. In their model, Thakur & Huang (2021) imposed (thus not driven by the model) a time-dependent shear modulus evolution to account for coseismic drop and postseismic recovery of elastic moduli. Erickson et al. (2017) has applied quasi-dynamic analysis to explore the effect of plasticity throughout the earthquake cycle. Preuss et al. (2020) have developed a 2.5D model with a visco-elasto-plastic crust subjected to rate- and state-dependent friction to model conjointly the rapid deformation in the elastic-brittle upper crust and the relaxation in the deeper viscoelastic crustal substrate and their influence on each other.

In all cases, irregardless of the end goal, setting the initial and boundary conditions (initial stresses for a single rupture, far field loading for seismic cycles) directly impacts both on-fault and off-fault dynamics. But what constraints do modellers have on the stress state of a fault? This will be discussed in the section below, followed-up by a review on how boundary conditions are set up in the community. We will then present a simple method to accurately define the initial stress field for in-plane conditions. To demonstrate its importance, we will explore its influence on-fault and off-fault deformation. Results are summarized and discussed in the last section.
# 2 Stress field and faulting

## 2.1 The Anderson theory

In the upper crust, large strains are accommodated by fault systems either seismically or aseismically. It either leads to the formation of new fractures in the crust, or reactivates previously existing faults. Fault systems evolve and acquire their general geometry by the progressive amalgamation of such fractures (Cowie & Scholz, 1992). In-situ measures of stresses on a fault, at any location, at any time are impossible to achieve, hence it has to be approached theoretically. A connection between the geometry of fault systems and the forces that formed them was first proposed by (Anderson, 1905). His theory related the initial formation of faults to the state of stress in the crust, under the assumption that rocks are isotropic, homogeneous and intact. This theory is based on the mathematical result that at every point of a stressed rock, three planes can be found on which no shear traction acts. Those planes are perpendicular to one another and are called the principal planes. The corresponding stresses acting along the three principal directions are called the principal stresses. By convention, the stresses are assumed to be positive in tension. Therefore, in a compressive regime, we have the following inequality:

$$ \sigma_1 < \sigma_2 < \sigma_3 \qquad (1) $$

Then Anderson assumed that, apart from those three planes, a plane with maximum tangential stress exists, on which the faulting should initiate. For brittle shear failure, and for a fluid saturated rock mass with pore fluid pressure $P_f$ (with $P_f > 0$), this Coulomb failure criterion may be written as:

$$ \tau_y = c + \mu(-\sigma_{eff}) = c + \mu(-\sigma_n - P_f) \qquad (2) $$

where $\sigma_n$ is the normal stresses acting on the plane and $\sigma_{eff} = \sigma_n + P_f$ corresponds to the effective normal stress. The variable $c$ is the cohesion and $\mu = \tan \phi$ is the coefficient of friction. The angle $\phi$ corresponds to the slope of the failure envelope on a Mohr diagram (Figure 1) and $\theta$, on the same diagram, corresponds to the angle between the normal to the fault plane and $\sigma_1$ (Figure 1). Hence we can write:

$$ \begin{aligned} \sigma_{yy} &= \left( \frac{\sigma_1 + \sigma_3}{2} \right) - \left( \frac{\sigma_3 - \sigma_1}{2} \right) \cos 2\theta \\ \sigma_{xx} &= \left( \frac{\sigma_1 + \sigma_3}{2} \right) + \left( \frac{\sigma_3 - \sigma_1}{2} \right) \cos 2\theta \\ \sigma_{xy} &= \left( \frac{\sigma_3 - \sigma_1}{2} \right) \sin 2\theta \end{aligned} \qquad (3) $$

The radius $R$ of the Mohr circle is then given by:

$$ R = \left( \frac{\sigma_3 - \sigma_1}{2} \right) = \sqrt{\left( \frac{\sigma_{xx} - \sigma_{yy}}{2} \right)^2 + \sigma_{xy}^2} \qquad (4) $$

In numerical studies, it is more common to define $\Psi$, the angle between the fault plane and $\sigma_1$ (Figure 1 & 2), such as $\Psi = \pi/2 - \theta$. Hence, as illustrated in Figure 1a, the relationship between the magnitudes of $\sigma_{yy}$ and $\sigma_{xx}$ depends on $\Psi$ in the following manner:

$$ \begin{aligned} |\sigma_{yy}| &\le |\sigma_{xx}| & \text{if } 0^\circ \le \Psi \le 45^\circ \\ |\sigma_{yy}| &> |\sigma_{xx}| & \text{if } 45^\circ < \Psi \le 90^\circ \end{aligned} \qquad (5) $$
a)

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

b)

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

c)

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

d)

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

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

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

<figure><img src="image_9.png" /><figcaption>Figure 1: a) Mohr–Coulomb criterion in (σ, τ)-space. b), c) & d) Orientation of the failure plane with respect to the largest principal stress for normal, strike-slip, and reverse faulting respectively. The angle φ correspond to the slope of the failure envelope on a Mohr diagram and defines the static friction μ = tan φ. The remaining angles Ψ, θ and ω are defined within the plane perpendicular to σ₂. Ψ and θ corresponds to the angle between σ₁ and the fault plane (Ψ) or, following the Mohr-Coulomb convention, the normal to the fault plane (θ). ω gives to the dip angle of the fault.</figcaption></figure>
Then, for an optimal angle $\theta_{opt} = \pi/4 + \phi/2$ shear failure occurs on a plane containing the $\sigma_2$ direction, when $\sigma_{yy} = \sigma_{eff}$ and $\sigma_{xy} = \tau_y$, i.e., when the failure envelope is tangential to the Mohr's circle. Relying on experimental studies, the internal friction $\mu$ generally lies between 0.5 and 1 (Jaeger and Cook, 1979). Hence the optimal angle $\theta_{opt}$ varies between 58° and 68° and the new fault should make an angle $\Psi$ with $\sigma_1$ between 32° and 22° respectively. Following the same reasoning for a pre-existing fault plane, if we consider a static coefficient of friction $\mu_s = 0.6$, which corresponds to a large variety of rocks and minerals (Byerlee, 1978), faults are optimally oriented when they make an angle $\Psi \simeq 30°$ with $\sigma_1$.

The foundation of the Anderson's theory lies in the observation that, due to Earth's free surface, it is imperative for one principal stress to be oriented subvertically, i.e., to be equal to $\sigma_v$. This gives rise to three fundamental stress regimes (Figure 1) depending on whether $\sigma_v = \sigma_1$, $\sigma_2$ or $\sigma_3$: reverse faulting ($\sigma_v = \sigma_3$), strike-slip faulting ($\sigma_v = \sigma_2$), and normal faulting ($\sigma_v = \sigma_1$). Keeping $\mu_s = 0.6$ as a reference value for the static friction, it is then expected to get sub-vertical strike-faults and a dip angle of $\omega \approx 30°$ and $\omega \approx 60°$ for a thrust and a normal fault respectively (Figure 1).

Despite the simplicity of the theory, observations are in accordance with the model for strike slip faults with low cumulative displacement (Anderson, 1951; Kelly et al., 1998). Earthquakes have also been registered on subvertical fault striking approximately 30° to the regional $\sigma_1$ with a subvertical $\sigma_2$: the 2000 Western Tottori earthquake in Japan (Sibson et al., 2012; Fukuyama et al., 2003; Yukutake et al., 2007) or the 2010 Mw 7.1 Darfield earthquake in New Zealand (Sibson et al., 2012). Moreover, borehole measurements and induced seismicity (Townend & Zoback, 2000), paleostresses inversion (Lisle et al., 2006) and earthquake focal mechanisms (Célérier, 2008) suggest that Andersonian state of stress prevails over large areas within the shallow crust. Of course exceptions to the theory exist too. Well-known examples, such as the San Andreas fault (Zoback et al., 1987), or the low-angle normal faults in Elba, Central Italy (Collettini & Sibson, 2001) or in the Cyclades Greece (Lecomte et al., 2010) are mis-oriented if we consider a coefficient of friction of 0.6 (Byerlee, 1978). But a lower friction on the fault plane (clay minerals) or a high pore pressure may well explain the discrepancy. If the fault essentially slips during earthquakes, weakening mechanisms such as the ones listed by Tullis & Schubert (2015) may well kick in (mineral breakdown, flash melting, thermal pressurization, etc...), leading to a very low effective coefficient of friction. Finally, polymodal faulting, with three or more sets of faults forming and slipping simultaneously are not compatible with the Anderson's assumption that faults form parallel to the intermediate principal stress, $\sigma_2$ (Healy et al., 2015). Hence, excluding these exceptions, it suggests that the Anderson's theory provides a useful framework for general considerations about fault generation and reactivation. And considering the lack of systematic, time-dependent in-situ stress measurements, probably the only one.

## 2.2 Non-exhaustive review on how initial stress state is prescribed in numerical models

In single dynamic rupture models, rupture grows under the control of the prescribed initial stresses. Theoretical analyses (Poliakov et al., 2002; Rice et al., 2005; Ngo et al., 2012) have illustrated the effect of initial stress field on the pattern of off fault damage activation (i.e., the potential failure area, and the orientation of secondary cracks) around a propagating crack. They show that the extent and location of secondary faulting (the activated zone) is strongly affected by the orientation of principal
stresses, set up by Ψ (the angle between σ₁ and the primary fault). Steep Ψ favors inelasticity on the extensional side and shallow Ψ on the compressional side. Moreover initial stresses seem to influence the potential for rupture to follow intersecting faults with different orientation rather than the primary fault (Kame et al., 2003b; Bhat et al., 2004; Fliss et al., 2005). Therefore, pre-stress orientation is a key parameter in numerical simulation of dynamic rupture with off-fault inelastic deformation.

Most of the numerical studies investigating the interactions between seismic rupture propagation and off fault damage usually set up a 2D planar strike-slip fault under plane-strain conditions (Andrews, 2005; Shi & Ben-Zion, 2006; Templeton & Rice, 2008; Dunham et al., 2011b; Thomas et al., 2017b; Okubo et al., 2019), such as displayed in Figure 1. Because of the 2D setting, the out of plane stress is often ignored when setting up the initial stress field. Usually, a normal stress and a shear stress are imposed, corresponding to two principal horizontal stresses σ₁ and σ₃ and an angle Ψ between fault and main stress direction. However, whether a fault is under a strike-slip stress field depends on the full 3D stress field (Figure 1). If the two in-plane stresses are σ₁ and σ₃, then the fault is in a proper strike-slip stress field. If the two in-plane stresses are σ₁ and σ₂, then the fault is in a reverse stress field. This would not change the slip on the fault because slip is restricted to the 2D plane. However, the ratio between the out of plane stress and the smallest principal in-plane stress would be smaller than expected under a proper 3D stress strike-slip stress field.

Despite the importance of initial stress field in modeling off fault deformation, the full three dimensional initial stress field and the importance of the out-of-plane stress, in 2D simulations, have never been discussed.

## 2.3 Criterion for accurate initial stress field in 2D plane-strain simulations

When considering a planar, strike-slip fault under plane-strain condition (Figure 2). The initial stress state, σᵢⱼ, (tensile positive) is represented as,

$$ \bar{\sigma} = \sigma_{ij} \equiv \begin{pmatrix} \sigma_{xx} & \sigma_{xy} & 0 \\ \sigma_{xy} & \sigma_{yy} & 0 \\ 0 & 0 & \sigma_{zz} \end{pmatrix} \qquad (6) $$

For convenience, we define the following ratios:

$$ \gamma = \frac{\sigma_{xx}}{\sigma_{yy}} \quad \& \quad \mu_0 = \frac{\sigma_{xy}}{-\sigma_{yy}} \qquad (7) $$

which leads to:

$$ \tan 2\Psi = -\tan 2\theta = -\frac{2\sigma_{xy}}{\sigma_{xx} - \sigma_{yy}} = \frac{2\mu_0}{\gamma - 1} \qquad (8) $$

Let σ₁, σ₂ and σ₃ represent the maximum, intermediate and minimum compressive principal stresses. For a strike-slip fault, σ₁ and σ₃ should lie on the x - y plane and σ₂ should be paral-
<figure><img src="image_10.png" /><figcaption>Figure 2: Modeling set up for a dynamic rupture on a strike-slip fault. a) Schematic of the fault zone in 3D. b) Zoom on the 2D plane that hosts the modeled rutpure. The orientation of maximum compressive stress $\sigma_1$ is set at an angle $\Psi$ to the fault plane that is governed by slip weakening friction law. The nucleation patch is set by increasing the initial shear stress $\sigma_{xy}$ slightly above the fault strength $\mu_s$ ($-\sigma_{yy}$) over a length $L_c$.</figcaption></figure>

Let to the z axis, i.e out of plane for a 2D simulation (Figure 1). Hence,

$$ 
\begin{align} 
\sigma_1 &= \left( \frac{\sigma_{xx} + \sigma_{yy}}{2} \right) - R = \left( \frac{\sigma_{xx} + \sigma_{yy}}{2} \right) \left( 1 + \sqrt{\left( \frac{\gamma - 1}{\gamma + 1} \right)^2 + \frac{4\mu_0^2}{(\gamma + 1)^2}} \right), \nonumber \\
\sigma_3 &= \left( \frac{\sigma_{xx} + \sigma_{yy}}{2} \right) + R = \left( \frac{\sigma_{xx} + \sigma_{yy}}{2} \right) \left( 1 - \sqrt{\left( \frac{\gamma - 1}{\gamma + 1} \right)^2 + \frac{4\mu_0^2}{(\gamma + 1)^2}} \right). \tag{9} \\
\sigma_2 &= \sigma_{zz} = \rho g z (1 - \lambda) \nonumber
\end{align} 
$$

where $\lambda$ is the pore pressure coefficient and $z$ is a parameter corresponding to the depth of the 2D slice where plane strain simulations are conducted. Under the plane-strain approximation, we also have the following relationship:

$$ \sigma_2 = \sigma_{zz} = \nu(\sigma_{xx} + \sigma_{yy}) = \nu\sigma_{yy}(\gamma + 1) \tag{10} $$

where $\nu$ is the Poisson's ratio.

We require that both $\sigma_1$ and $\sigma_3$ should be compressive, i. e. $\sigma_1 < \sigma_3 < 0$. This implies that

$$ \sqrt{\left(\frac{\gamma - 1}{\gamma + 1}\right)^2 + \frac{4\mu_0^2}{(\gamma + 1)^2}} < 1 \tag{11} $$

since $\sigma_{xx} + \sigma_{yy} < 0$. The principal stress state must also satisfy the inequality given in equation 1.

$$ 
\begin{align} 
\sigma_1 < \sigma_2 &\Rightarrow \left( \frac{\sigma_{xx} + \sigma_{yy}}{2} \right) \left( 1 + \sqrt{\left( \frac{\gamma - 1}{\gamma + 1} \right)^2 + \frac{4\mu_0^2}{(\gamma + 1)^2}} \right) < \nu(\sigma_{xx} + \sigma_{yy}) \nonumber \\
&\Rightarrow \sqrt{\left( \frac{\gamma - 1}{\gamma + 1} \right)^2 + \frac{4\mu_0^2}{(\gamma + 1)^2}} > 2\nu - 1 \tag{12}
\end{align} 
$$
This is trivially satisfied as $\nu$ is always smaller than 0.5.

$$ \sigma_3 > \sigma_2 \Rightarrow \left(\frac{\sigma_{xx} + \sigma_{yy}}{2}\right) \left(1 - \sqrt{\left(\frac{\gamma - 1}{\gamma + 1}\right)^2 + \frac{4\mu_0^2}{(\gamma + 1)^2}}\right) > \nu(\sigma_{xx} + \sigma_{yy}) \quad (13) \\ \Rightarrow \sqrt{\left(\frac{\gamma - 1}{\gamma + 1}\right)^2 + \frac{4\mu_0^2}{(\gamma + 1)^2}} > 1 - 2\nu $$

Thus, the stress field for a strike-slip fault must satisfy the following criterion:

$$ 1 - 2\nu < \sqrt{\left(\frac{\gamma - 1}{\gamma + 1}\right)^2 + \frac{4\mu_0^2}{(\gamma + 1)^2}} < 1 \quad (14) $$

This inequality is plotted in Figure 3 as a function of $\gamma$ and $\mu_0$. The shaded areas define what we will from now on refer to as the “correct” regime (strike-slip stress field). The white areas represent (1) the “forbidden” regime (reverse faulting stress field) for which, this criterion is violated and (2) the stress field corresponding to a tensile regime. Superimposed on this graph are the initial parameters for several studies modeling strike-slip motion using plane-strain approximation. Three out of six studies used a far field stress field that favors reverse faulting, i.e the out-of-plane stress is $\sigma_3$ and not $\sigma_2$ as it should be (Figure 3a). If the modelling is only performed on a 2D plane, like the vast majority of the published studies, the fault will still have a strike-slip motion even if this condition is not satisfied simply because the motion is restricted along one plane. However, the whole stress field would favor reverse faulting, which will impact any model of off-fault deformation as we will demonstrate soon.

In Figure 3c, we plot the same criterion as a function of $\Psi$ and the seismic ratio $S$, as defined by (Andrews, 1976; Das & Aki, 1977),

$$ S = \frac{\mu_s(-\sigma_{yy}) - \sigma_{xy}}{\sigma_{xy} - \mu_d(-\sigma_{yy})} = \frac{\mu_s - \mu_0}{\mu_0 - \mu_d} \quad (15) $$

where $\mu_s$ and $\mu_d$ correspond to the static and dynamic coefficient of friction respectively, for a slip weakening law. $S$ helps determine whether a rupture in 2-D, under homogeneous conditions, can reach supershear velocities ($S < 1.77$), or remains sub-Rayleigh ($S > 1.77$). For a fixed value of $\Psi$, fulfilling the criterion for strike-slip faulting strongly depends on $\mu_d$. Interestingly, setting up a proper strike-slip stress field in agreement with Anderson theory, i.e. $\Psi \simeq 30^{\circ}$ for $\mu_s = 0.6$, requires supershear parameters, and this holds true for a large range of realistic $\mu_d$ values.

# 3 Methods

## 3.1 Numerical model setup

In this study, we explore the boundary between the forbidden regime and the correct regime by modeling rupture on a 1D right-lateral fault in a 2D medium under plane-strain approximation (Figure 2). We particularly focus on the influence of the initial stress conditions on off-fault stresses and inelastic deformation. We use the 2-D spectral element code SEM2DPACK (Ampuero, 2012). Rupture propagation along the fault plane is governed by a slip-weakening friction law (e.g., Palmer & Rice, 1973). Slip occurs when the on-fault shear stress reaches the shear strength $\tau_f = \mu^*(-\sigma_{yy})$. The
<figure><img src="image_11.png" /><figcaption>Figure 3: Criterion for accurate initial stress field under plane-strain approximation displayed as $\gamma - \mu_0$, plots (a) & (b) and $\Psi - S$ for plot (c). The shaded areas represent the conditions for which the initial stresses favor a strike-slip motion, i.e. when equation (14) is satisfied. In the $\gamma - \mu_0$ plot, this criterion depends on $\nu$ (0.25 and 0.27 for plot (a) and (b) respectively), but is independent from other parameters. Contours lines shows $\Psi$ values. Color dots shows the initial parameters for published simulations of strike-slip faulting with plane-strain approximation. In plot (b), the shaded red areas shows the two sets of parameters explored in this study, with $\Psi = 15^\circ$ and $\Psi = 55^\circ$. The red empty circles correspond to the end members. In subplot (c) the criterion is represented as a function of $S$ and $\Psi$ for $\mu_s = 0.6$. Plotting the area corresponding to an accurate initial stress depends on a third parameter, here given by $\mu_d$. We plot the criterion for two cases, $\mu_d = 0.09$ (dotted area) and $\mu_d = 0.34$ (gray area).</figcaption></figure>
<table><caption>Table 1: Constants used in all models</caption><thead><tr><td>Parameter</td><td>Symbol</td><td>Value</td></tr></thead><tbody><tr><td>depth</td><td>z</td><td>-2.5 km</td></tr><tr><td>material density</td><td>ρ</td><td>2700 kg.m<sup>-3</sup></td></tr><tr><td>S-wave velocity</td><td>c<sub>s</sub></td><td>3120 km.s<sup>-1</sup></td></tr><tr><td>P-wave velocity</td><td>c<sub>p</sub></td><td>5600 km.s<sup>-1</sup></td></tr><tr><td>Poisson's ratio</td><td>ν</td><td>0.27</td></tr><tr><td>characteristic slip</td><td>D<sub>c</sub></td><td>1 m</td></tr><tr><td>pore pressure coefficient</td><td>λ</td><td>0.4</td></tr></tbody></table>

friction coefficient $\mu^*$ depends on the cumulated slip ($\delta$) and drops from a static $\mu_s$ to a dynamic $\mu_d$ value over a characteristic distance ($D_c$). The rupture is artificially nucleated within a patch where the initial shear stress $\sigma_{xy}$ is set just above the fault strength (Figure 2). Further on, since the difference in nucleation duration $t_{nuc}$ between two models has no physical meaning, we will shift all the results in time so $t_{nuc} = t_0$. Following Kame et al. (2003a), the minimum nucleation size $L_c$ determined by the energy balance for a slip weakening law is:

$$L_c = \frac{16}{3\pi} \frac{\mu G}{(\sigma_{xy} - \tau_r)^2} \qquad (16)$$

where G is the fracture energy, defined as:

$$G = \frac{1}{2} D_c (\tau_p - \tau_r) \qquad (17)$$

and where $\tau_p = \mu_s(-\sigma_{eff})$ and $\tau_r = \mu_d(-\sigma_{eff})$ are respectively the peak and residual stresses.

We further use the process zone size $R_0$ for a quasi-stationary crack, to normalize the length scales in our results. It corresponds to the length over which the friction drops, with ongoing slip, from the peak strength to the residual strength. Following Day et al. (2005) it is given by,

$$R_0 = \frac{9\pi}{32(1-\nu)} \frac{D_c \mu}{(\mu_s - \mu_d)(-\sigma_{yy})} \qquad (18)$$

In our models, $R_0$ is approximately equal to 1.0 km for $\Psi = 55^\circ$ and 1.6 km for $\Psi = 15^\circ$ (see section 3.2 for an explanation on the parameters). We use a resolution of 30 m to ensure that the problem is correctly resolved. In order to scale the problem, all the modelled faults have a length of $30R_0$.

## 3.2 Initial parameters

In order to compare our study to the available literature, we run two sets of models with $\Psi = 55^\circ$ and $\Psi = 15^\circ$ respectively (Figure 3a). To explore the effect of initial stresses on off-fault and on-fault deformation, for each set of models, we run three simulations with the correct strike-slip set-up (cases i, ii, iii for $\Psi = 55^\circ$ and cases vii, viii, ix for $\Psi = 15^\circ$) and three within the so-called forbidden regime (cases iv, v, vi for $\Psi = 55^\circ$ and cases x, xi, xii for $\Psi = 15^\circ$).

To achieve that goal, we adopt the following strategy. Some parameters are kept constant between
all the simulations: the elastic properties, the depth, the characteristic slip in the friction law and a hydrostatic pore pressure condition (Table 1). We have two $\Psi$ that define the orientation of maximum compressive stress $\sigma_1$ with respect to the fault. We set the S ratio (equation 15), favouring a value that would lead to subshear rupture ($S = 2$) for $\Psi = 55^{\circ}$. However, for $\Psi = 15^{\circ}$, this lead to a small dynamic stress drop which prevents the rupture from propagating (dying cracks). Hence we choose to set the ratio to $S = 1$ for $\Psi = 15^{\circ}$ (the rupture can evolve to a supershear earthquake). We further fix the stress drop ($\sim 10$ MPa):

$$ \Delta\tau = (\mu_0 - \mu_d)(-\sigma_{yy}) \qquad (19) $$

and, in agreement with laboratory values (Jaeger, 1979), we make $\mu_s$ varies between 0.52 and 0.68 to get the six models defined above.

The others parameters can be determined using the following set of equations. In hydrostatic pore pressure condition, the vertical, out-of-plane stress $\sigma_{zz}$ is given by

$$ \sigma_{zz} = \rho g z (1 - \lambda) \qquad (20) $$

We then need to compute the ratio $\gamma$:

$$ \begin{array}{l} \text{eq. 7 & 10} \quad \Rightarrow \gamma = \frac{\sigma_{zz}}{\nu\sigma_{yy}} - 1 \\ \text{eq. 19 & 15} \quad \Rightarrow \gamma = \frac{\sigma_{zz}(\mu_0 - \mu_s)}{\nu S \Delta\tau} - 1 \\ \text{eq. 8} \quad \Rightarrow \gamma = \frac{\sigma_{zz}(0.5(\gamma - 1)\tan 2\Psi - \mu_s)}{\nu S \Delta\tau} - 1 \qquad (21) \\ \Rightarrow \gamma = \frac{(\alpha + \beta) + \mu_s\sigma_{zz}}{(\beta - \alpha)} \\ \text{with } \alpha = \nu S \Delta\tau \ \& \ \beta = 0.5\sigma_{zz} \tan 2\Psi \end{array} $$

Knowing $\gamma$, we can derive:

$$ \mu_0 = \frac{(\gamma - 1)}{2} \tan 2\Psi \qquad (22) $$

$$ \mu_d = \mu_0 - \frac{\mu_s - \mu_0}{S} \qquad (23) $$

$$ \sigma_{yy} = \frac{\sigma_{zz}}{\nu(\gamma - 1)} \qquad (24) $$

$$ \sigma_{xx} = \gamma\sigma_{yy} \qquad (25) $$

$$ \sigma_{xy} = \mu_0(-\sigma_{yy}) \qquad (26) $$

The values of initial parameters for the end-member models are summarized in Table 2 while the corresponding Mohr-Coulomb diagrams illustrating their initial stress fields can be found in the supplementary materials (see Figure S1). For a fixed angle $\Psi$, note that because we set S constant between the six models, a constant stress drop is equivalent to a constant fracture energy (equation 17). Therefore, the characteristic length scales for the friction law, $R_0$ and $L_c$, are also constant.
<table><caption>Table 2: Initial parameters used for the two end-members of each set of models.</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Symbol</th><th colspan="4">Initial stress field</th></tr><tr><th>case (i)<br>strike-slip</th><th>case (vi)<br>reverse</th><th>case (vii)<br>strike-slip</th><th>case (xii)<br>reverse</th></tr></thead><tbody><tr><td rowspan="4">Input</td><td>angle</td><td>Ψ</td><td>55<sup>o</sup></td><td>55<sup>o</sup></td><td>15<sup>o</sup></td><td>15<sup>o</sup></td></tr><tr><td>seismic ratio</td><td>S</td><td>2</td><td>2</td><td>1</td><td>1</td></tr><tr><td>stress drop (MPa)</td><td>Δτ</td><td>10.2</td><td>10.2</td><td>10.6</td><td>10.6</td></tr><tr><td>static friction</td><td>μ<sub>s</sub></td><td>0.68</td><td>0.52</td><td>0.68</td><td>0.52</td></tr><tr><td rowspan="12">Resulting parameters</td><td>dynamic friction</td><td>μ<sub>d</sub></td><td>0.34</td><td>0.16</td><td>0.19</td><td>0.09</td></tr><tr><td rowspan="4">full stress tensor (MPa)</td><td>σ<sub>yy</sub></td><td>-91</td><td>-84</td><td>-43</td><td>-50</td></tr><tr><td>σ<sub>xx</sub></td><td>-61</td><td>-67</td><td>-108</td><td>-102</td></tr><tr><td>σ<sub>zz</sub></td><td>-42</td><td>-42</td><td>-42</td><td>-42</td></tr><tr><td>σ<sub>xy</sub></td><td>41</td><td>23</td><td>19</td><td>15</td></tr><tr><td rowspan="3">principal stresses (MPa)</td><td>σ<sub>1</sub></td><td>-120</td><td>-101</td><td>-113</td><td>-106</td></tr><tr><td>σ<sub>2</sub></td><td>-42</td><td>-51</td><td>-42</td><td>-45</td></tr><tr><td>σ<sub>3</sub></td><td>-32</td><td>-42</td><td>-38</td><td>-42</td></tr><tr><td>process zone (m)</td><td>R<sub>0</sub></td><td>1047</td><td>1047</td><td>1570</td><td>1570</td></tr><tr><td>nucleation length (m)</td><td>L<sub>c</sub></td><td>6562</td><td>6562</td><td>4374</td><td>4374</td></tr></tbody></table>

# 4 Results

## 4.1 Role of pre-stresses in determining the yield criterion

In this section, we run simulations with an elastic medium and we compute different yield criterion commonly used in the literature to determine the inelastic deformation. It is important to note that in those cases, the medium has a pure elastic behavior and the off-fault plastic deformation is calculated a posteriori.

### 4.1.1 Role of the out-of-plane stress in computing the plastic yield criterion
We first examine the importance of accounting for the out-of-plane stress in determining the off-fault plastic deformation. Here we compute the Drucker-Prager criterion:

$$F_{DP} = \sqrt{J_2} + \mu_s p \qquad (27)$$

where $p \equiv I_1/3 = (\sigma_{xx} + \sigma_{yy} + \sigma_{zz})/3$ is the hydrostatic stress derived from the first invariant $I_1$ of the stress tensor and $J_2 = s_{ij}s_{ij}/2$ corresponds to the second invariant of the deviatoric stress tensor (with $s_{ij} = \sigma_{ij} - p\delta_{ij}$). However, in many of the published 2D studies (Templeton & Rice, 2008; Dunham et al., 2011b, among others) the out-of-plane stress, $\sigma_{zz}$, is assumed to be the mean of the in-plane stresses i.e. $\sigma_{zz} = (\sigma_{xx} + \sigma_{yy})/2$. This particular choice of $\sigma_{zz}$ makes the Mohr-Coulomb and Drucker-Prager yield surfaces coincide in 2D. Therefore, the invariants can be computed using the in-plane stress tensor components i. e. $p^{ip} = (\sigma_{xx} + \sigma_{yy})/2$ which consequently changes the second invariant of the deviatoric stress tensor as well (further referred as $J_2^{ip}$). We note the Drucker-Prager criterion, solely using the in-plane stresses as follow:

$$F_{DP}^{ip} = \sqrt{J_2^{ip}} + \mu_s p^{ip} \qquad (28)$$
To illustrate the contribution of the out-of-plane stress in computing the invariants we use case (i) with $\Psi = 55^\circ$ and $S = 2$ (see Table 2). We compare the hydrostatic stress and the square root of the second invariant of the deviatoric stress tensor obtained with and without the out-of-plane stress (Figures 4a & b) when the fault has ruptured about 5 times the process zone $R_0$. We can observe that in both cases, the hydrostatic stress is higher in the compressional quadrants. However, if the out of plane stress is ignored, the hydrostatic stress is overestimated (stresses are positive in tension) by up to 20 MPa (or two times the dynamic stress drop) as displayed in Figure 4c. On the other hand, the second invariant of the deviatoric stress tensor is underestimated if we ignore $\sigma_{zz}$ (Figures 4d, e & f). The difference between $J_2$ and $J_2^{ip}$ is up to ~1.5 times the dynamic stress drop. This results in a significant difference in the estimation of the Drucker-Prager criterion (Figure 5a, b & c). In both cases, $F_{DP}$ or $F_{DP}^{ip}$ are positive in the tensional quadrants. However, taking into account the out-of-plane stress not only increases the area where the Drucker-Prager criterion is positive, i.e., where the plastic deformation is expected, but the overall magnitude of the plastic strain is also higher. Hence this leads to an underestimation of the plastic deformation. As expected, Figure S2 reveals consistent findings when computing the Drucker-Prager criterion, solely using the dynamic change of stresses due to the rupture propagation.

Changing the initial state of stress, using case (vii), i.e. for $\Psi = 15^\circ$ and $S = 1$ (see Table 2) we obtained even more drastic differences. When the out-of-plane stress is ignored (Figure 5e) $F_{DP}^{ip}$ is pretty much negative everywhere, which may lead to the interpretation that no plastic deformation is happening. Whereas, when accounting for $\sigma_{zz}$, even if the magnitude is about four times smaller than for case (i), we record positive $F_{DP}$, notably within the compressional quadrants.

### 4.1.2 Role of the reverse versus strike-slip stress field on the plastic yield criterion
Now that we have emphasized the importance of including the full stress tensor in computing the Drucker-Prager criterion, we look at the influence of setting up a proper strike-slip initial stress field, ensuring that the pre-stress conditions satisfy $\sigma_1 < \sigma_2 = \sigma_{zz} < \sigma_3$.

In Figure 6, for $\Psi = 55^\circ$, $S = 2$ and for $\Psi = 15^\circ$, $S = 1$ we compare the Drucker-Prager criterion of the end-member models, i.e case (i), with case (vi) and case (vii) with case (xii). The pre-stress conditions for cases (i) and (vii) correspond to a strike-slip fault, whereas cases (vi) and (xii) have the initial stress field of a reverse fault. The yield criterion is computed when the fault has ruptured about 5 times the process zone $R_0$. Note that, due to the difference in S ratio for the simulations with $\Psi = 55^\circ$ and $\Psi = 15^\circ$, it is not meaningful to compare the models with different angles, as the dynamic rupture is very different, and by extent, the off-fault deformation.

For initial reverse stress conditions, when $\Psi = 55^\circ$ (Figure 6b, case vi), positive $F_{DP}$ is observed in the tensile quadrants. When $\Psi = 15^\circ$ (Figure 6e, case xii), positive $F_{DP}$ is observed in the compressive quadrants, with a magnitude lower than for case (vi). When the initial stress field is properly set-up (Figure 6a & d), the areas recording a positive $F_{DP}$ are much larger and the absolute value is also higher, as illustrated by Figure 6c & f. Since the models are run within a purely elastic medium and since the rupture dynamics are very similar, the differences are essentially linked to the background stresses.

Figure 7 shows the continuum of all models for $\Psi = 55^\circ$. Unlike Figure 6 for which we plot $F_{DP}$ at one particular time step, here we plot the maximum value of $F_{DP}$ induced in the off-fault medium by
<figure><img src="image_12.png" /><figcaption>Figure 4: First invariant of the stress tensor (a,b) and second invariant of the deviatoric stress tensor (d,e) computed with the 3D stress tensor (a,d), the in-plane stress tensor (b,e), at t = 2.9 seconds for case (i) with S = 2 and Ψ = 55 °. Figures (c,f) give respectively the difference between the way of computing the first invariant of the stress tensor and the second invariant of the deviatoric stress tensor. Invariants are normalized by the dynamic stress drop Δτ (equation 19). Slip rate on the fault (black curves) is super-imposed on the graphs.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 5: Drucker-Prager criterion for case (i) with $S = 2$ and $\Psi = 55^\circ$ (a,b) and for case (vii) with $S = 1$ and $\Psi = 15^\circ$ (d,e), at $t = 2.2$ seconds. We compute the invariants accounting for $\sigma_{zz}$ (a,d), or only using the in-plane stress field (b,e). We show the differences $F_{DP} - F_{DP}^{ip}$ in Figures (c,f) for cases (i) and (vii), respectively. The difference is computed for the positive values only. Different versions of the Drucker-Prager criteria and their differences are normalized by the dynamic stress drop $\Delta\tau$ (equation 19). Slip rate on the fault (black curves) is super-imposed on the graphs.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 6: Drucker-Prager criterion for cases (i) & (vi), with $S = 2$ and $\Psi = 55^\circ$ (a,b) and for cases (vii) & (xii), with $S = 1$ and $\Psi = 15^\circ$ (d,e). Cases (i) & (vii) correspond to an initial strike-slip stress field (a,d) and cases (vi) & (xii) to an initial reverse stress field (b,e). Figures (c,f) give the difference between cases (i) & (vi), and between cases (vii) & (xii), respectively. The four models are plotted against the criterion for accurate initial stress in Figure 3. Drucker-Prager criterions and differences are normalized by the dynamic stress drop $\Delta\tau$ (equation 19). Slip rate on the fault (black curves) is super-imposed on the graphs.</figcaption></figure>

the full rupture. Crossing the boundary between initial strike-slip stress field and reverse pre-stress conditions does not change the results qualitatively. However, the area with positive $F_{DP}$ is larger in the strike-slip case, and we observe higher absolute values, with up to twice the dynamic stress drop.

### 4.1.3 Role of the initial stress field for off-fault rupture modes
We have shown that the initial stress conditions influence the values of $F_{DP}$ criterion induced by the rupture. Yet the Drucker-Prager criterion uses the invariants of the stress tensor, which does not predict the preferential orientation for failure, a useful information in seismic risk assessment for example. Another criterion often used in the literature is therefore the Coulomb stress change due to the rupture on the main fault (Stein et al., 1997; King et al., 1994; Thomas et al., 2017a; Canitano et al., 2021).

$$ \Delta CFF = (\Delta\tau_y + \mu_s\Delta\sigma_{eff}) \qquad (29) $$

where $\Delta\sigma_t$ and $\Delta\sigma_n$ are the change of shear and normal stress respectively on the optimally oriented plane, induced by the seismic rupture for a given direction $\theta$ with respect to the principal stress $\sigma_1$ (in relation to the Mohr-Coulomb circle displayed in Figure 1). In Figure 8 we explore the same cases
<figure><img src="image_15.png" /><figcaption>Figure 7: Drucker-Prager criterion for all models with $\Psi = 55^\circ$ and $S=2$. We plot the maximum value $F_{DP}$ between $t=0$ and $t=4$ seconds. Only positive values are considered (hence the color discontinuity.) We display the results from the strike-slip end-member (upper left corner) to the reverse end-member (lower right corner). The initial parameters are shown in the upper right inset, in the same representation as in Figure 3a. The lower left inset shows the difference between the two end members. Drucker-Prager criterions and differences are normalized by the dynamic stress drop $\Delta\tau$ (equation 19). Slip rate on the fault (black curves) is super-imposed on the graphs.</figcaption></figure>
that in Figure 6, i.e., the two end-members for $\Psi = 55^{\circ}$, S = 2 (cases i & vi) and for $\Psi = 15^{\circ}$, S = 1 (cases vii & xii). The background color corresponds to local values of $\Delta CFF$. The yield criterion is computed when the fault has ruptured about 5 times the process zone $R_0$. As observed previously with the Drucker-Prager criterion, areas with positive Coulomb stress change are larger and with higher $\Delta CFF$ values for $\Psi = 55^{\circ}$ Figures 8a & b. Likewise, the areas likely to induce off-fault deformation are larger when the initial stress-field is properly set up, and this for the two tested values of $\Psi$ Figures 8a & c.

On top of these snapshots, we compute the local preferential orientations for failure and the corresponding type of faulting induced. An off-fault strike-slip failure means that locally, $\sigma_1$ and $\sigma_3$ are in-plane. If $\sigma_3$ is out-of-plane, the local preferred type of failure will be that of a reverse fault. Figure 8a shows that if the initial strike-slip stress field is correctly set up, the rupture induces only strike-slip off-fault failures. However, if the initial stress field actually corresponds to that of a reverse faulting, the outcomes are different (Figure 8b). The main rupture (strike-slip by default since we are in 2D) only influences the stress field close to the rupture tips. Far from the fault and within the nucleation area, the propagating rupture does not control the local stress field, but its initial value does. Hence we observe off-fault reverse failures for case (vi).

## 4.2 Dynamic simulation with inelastic rheology

We have shown, using linear elasticity, that the pre-stress conditions can significantly affect the assessment of the different yield criteria used to estimate the off-fault deformation. In this section, we now investigate the role of the initial stress field on the dynamically triggered off-fault deformation, and its feedback on the seismic rupture, using different modelling strategies.

### 4.2.1 Role of the initial stress field in plastic deformation
We first use the off-fault plasticity model implemented in SEM2DPack following Andrews (2005). The inelastic response of the medium is characterised by a Coulomb criterion using in-plane stresses with $\mu = 0.75$ and $c = 30$ MPa so that the initial stress state of the medium is below the yield criterion. The maximum shear stress over all orientations is:

$$ \tau = \sqrt{\sigma_{xy}^2 + [(\sigma_{xx} - \sigma_{yy})/2]^2} \qquad (30) $$

The Coulomb criterion is:

$$ \tau \le c \cos \mu - 0.5(\sigma_{xx} + \sigma_{yy}) \cos \mu \qquad (31) $$

At each time step of the calculation, stress components are first incremented elastically. Then, if the Coulomb criterion is violated, stress components are recomputed so that a part of the deformation is accommodated inelastically (see Andrews (2005); Duan & Day (2008) for details about the method).

We compute the plastic deformation induced by the rupture for the two end-members cases (i) and (vi) previously studied ($\Psi = 55^{\circ}$ and S=2). In Figure (9), the results are illustrated by the second invariant of the deviatoric strain tensor $\gamma_t$:

$$ \gamma_t = \sqrt{2e_{ij}e_{ij}} \qquad (32) $$
<figure><img src="image_16.png" /><figcaption>Figure 8: Coulomb stress change ΔCFF due to the dynamic rupture, computed on optimally oriented planes and normalized by the dynamic stress drop Δτ (equation 19), for cases (i) & (vi), with S = 2 and Ψ = 55° (a,b) and for cases (vii) & (xii), with S = 1 and Ψ = 15° (c,d). Cases (i) & (vii) correspond to an initial strike-slip stress field (a,c) and cases (vi) & (xii) to an initial reverse stress field (b,d). Superimposed on this snapshot are the conjugates planes that give a maximum value of ΔCFF. The symbols inform on the expected mode of rupture: two crossed lines represent the two conjugate planes for strike-slip ruptures, the beach ball with a dark center gives the orientation of the conjugate planes for reverse faulting, following the classic seismological convention. Those symbols are displayed in the areas where ΔCFF > 0.2. Slip rate on the fault (grey curve) is super-imposed on the graphs.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 9: Plastic deformation for the models with $\Psi = 55^{\circ}$, S=2, with initial strike-slip stress field for case (i) and reverse stress field for case (vi). Here we show the second invariant of the deviatoric strain tensor $\gamma$. Slip rate on the fault (black curve) is super-imposed on the graphs.</figcaption></figure>

where $e_{ij} = \epsilon_{ij} - \frac{1}{3}\epsilon_{kk}\delta_{ij}$. The pattern of off-fault deformation is significantly different between the two simulations. The extent of the plastic deformation is much larger for an initial strike-slip stress state (Figure 9a), as expected based on the results of section 4.1. This induces differences in the rupture dynamics with a decrease of rupture speed and slip rate (Figure 10d) and a lower cumulative slip (Figure 10c), compared to case (vi) that has an initial stress field that corresponds to reverse faulting. Therefore, when off-fault inelastic deformation is taken into account, even modest differences in the initial stress condition affect significantly both the evolution of the off-fault medium and the slip dynamics.

We note that in a non-cohesive medium ($c=0$), those differences are even more emphasised (Figure S3a & b). For initial strike-slip stress conditions, the rupture decelerates rapidly, its propagation is prevented by the intense inelastic deformation of the medium. The off-fault deformation is localised and optimally oriented with respect to the far field stress orientation. It is equivalent to the creation of a new optimally oriented fault. Hence, it is interesting to note that modelling off-fault deformation for a 2-D strike-slip fault, with the appropriate initial stress field, requires a certain cohesion in order to fully rupture the prescribed fault plane. We also note that for a non-cohesive medium, when $\Psi = 15^{\circ}$ and S=1 (Figure S3c & d), little deformation occurs for the strike-slip case (vii) and none for the reverse case (xii). When $\Psi = 30^{\circ}$ (optimally oriented fault, Figure S3e & f), deformation only occurs on the main fault. Therefore, a flat fault has to be mis-oriented to produce significant off-fault deformation.
<figure><img src="image_18.png" /><figcaption>Figure 10: Cumulative slip (a-c-e-g) and slip rate (b-d-f-h) plotted every 0.4 seconds for $\Psi = 55^\circ$ and $S = 2$, with initial reverse stress fields (dark solid lines, case *vi*) and initial strike-slip stress field (dotted color lines, case *i*).</figcaption></figure>
### 4.2.2 Role of the initial stress field on dynamic damage
The last set of simulations use a micromechanics-based model to determine the dynamically-triggered off-fault damage and its feedback on the rupture dynamic. Inelastic deformation can occur in the model by either crack opening or crack propagation from initial flaws. Using an energy-based approach, at each time step, the corresponding change in elastic moduli, and hence the constitutive law, is determined (see Thomas et al. (2017b) for constitutive equations). The current inelastic state of the medium is characterized by the scalar D, the fraction of volume occupied by microcracks:

$$D = \frac{4\pi}{3} N_v (a \cos \Phi_c + l)^3 \quad (33)$$

where a is the initial microcrack radius, l, the wings cracks length as they are growing parallel to $\sigma_1$ ($l = 0$ at $t = 0$) (see figure 11) and $N_v$ the volume density of cracks. Initial flaws are all aligned at the same optimal angle $\Phi_c = \frac{1}{2} \tan^{-1}(1/\mu_c)$ to $\sigma_1$, with $\mu_c$ being the friction coefficient for the microcracks. In their model, Bhat et al. (2012b) derive the crack growth law by comparing the stress intensity factor at the microcrack tips to the experimentally-determined initiation and propagation toughness.

$D$ varies between 0 and 1, the maximum value corresponding to the coalescence stage that leads to the macroscopic fracture of the solid. See Bhat et al. (2012b) and Thomas et al. (2017b) for further details on the method.

Figure 11 shows the damage density induced by the rupture for the end-member cases (i) and (vi) with $\Psi = 55^\circ$ and S=2. The distribution of pre-existing flaws is homogeneous ($D = 0.1$ at $t = 0$) in Figures 11a & b. For the models in Figures 11c & d we assume an exponential decay of initial damage with fault normal distance, as described in several field studies (e.g., Vermilye & Scholz, 1998; Wilson et al., 2003; Mitchell & Faulkner, 2009). The initial damage density varies from $D = 0.5$ to $D = 0.1$ over a distance equivalent to the process zone $R_0 \sim 1$ km. The characteristic scale was selected to match the width of the damage zone that may have been created by past earthquakes, in this case defined in Figure 11a. In all scenarios, to prevent off-fault damage at the beginning of the simulation due to far field loading, we set the friction on the microcracks, $\mu_c$ to 0.75. Hence the observed damage is dynamically triggered by the seismic rupture.

Similar to the models discussed earlier, for these particular stress states, damage essentially occurs in the tensional quadrants. The damage zone is also significantly wider for initial strike-slip stress conditions (Figures 11a & c). Previous studies, in comparison to simulations with a pure elastic medium, have underlined the effect of damage on slip rate and rupture velocity (slow down) and to a lesser extent the cumulative slip (Thomas et al., 2017b; Thomas & Bhat, 2018). Here, when comparing in Figure 10 the models with initial reverse stress fields (dark solid lines) and strike-slip stress field (dotted color lines), we observe differences in cumulative slip (e & g), slip rate and rupture velocity (f & h), both for models with an uniform medium or an initial damage zone. That is because, unlike for the reverse case, when the initial stress field is correctly set up, damage occurs ahead of the rupture tip, thus changing the P- and S-wave speeds in the medium, which ultimately slows down the rupture velocity. This effect is even more pronounced when the earthquake ruptures a fault with a pre-existing damage zone (Figures 10g & h).
<figure><img src="image_19.png" /><figcaption>Figure 11: Simulation of a dynamic rupture with off-fault damage within an homogeneous medium (a,b) and within a medium a pre-existing damage zone, with an exponential decay of damage density away from the fault (c,d). We explore the two end-member cases for $S = 2$ and $Ψ = 55°$: case (i) corresponds to an initial strike-slip stress field (a,c) and cases (vi) to an initial reverse stress field (b,d). The colors represent the density of microcracks in the medium. Slip rate on the fault (red curve) is super-imposed on the graphs. The inset shows a schematic representation of the initial crack distribution and of the crack geometry (after Thomas et al. (2017b)).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 12: Synthesis of the effect of initial stress on damage zone width for cases (i) and (vi) with $\Psi = 55^\circ$ and S=2, when the rupture has reached seven times the static process zone size. Results from the different models are displayed in red for case (vi), with reverse pre-stress conditions, and in green for case (i), with strike-slip pre-stress conditions. a) Expected damage zone width using purely elastic models. We use the Coulomb stress change with a threshold of $\Delta CFF = 0.3\Delta\tau$. b) Damage zone width computed with the model following Andrews (2005), using a threshold value of plastic deformation of 0.01%. c) Damage zone width computed with the micromechanical model, using a threshold value of D = 0.3 (medium with a uniform initial damage density $D_0 = 0.1$). The shaded areas correspond to $\pm 10\%$ changes in the threshold value.</figcaption></figure>

# 5 Discussion and Conclusions

As discussed above, in the last two decades, the wealth of observations have underlined the importance of off-fault deformation and complex structure in fault zone behavior. Numerical models have been developed to incorporate these key ingredients (Andrews, 2005; Shi & Ben-Zion, 2006; Templeton & Rice, 2008; Dunham et al., 2011b; Thomas et al., 2017b; Okubo et al., 2019, among others). However, due to numerical limitations, they have been developed mostly in 2D. As a consequence, when setting up the initial stress field, only the in-plane stresses are defined and the out of plane stress is often ignored, or assumed to be the mean of the in-plane stresses.

In several studies, the effect of initial stress state have been investigated in terms of the orientation of $\sigma_1$ with respect to the fault direction. They illustrated the significant influence of $\Psi$ on the pattern of off-fault damage (Poliakov et al., 2002; Rice et al., 2005; Ngo et al., 2012; Templeton & Rice, 2008). However, the impact of the relative importance of the principal stresses, at constant angle $\Psi$, has not been discussed. In this study, we run 2D plane-strain simulation of a strike slip faulting to illustrate the key role of the 3D faulting regime on off-fault deformation.

Using an elastic medium, we first show that, even if the initial stress field is rightfully set up, ignoring the out-of-plane stress (here $\sigma_{zz}$) leads to an underestimation of the inelastic deformation, both in extend and magnitude. If the initial stress field is on top wrongly defined (reverse faulting), ignoring $\sigma_{zz}$ will on contrary lead to an over-estimation of the inelastic deformation. Then, we demonstrate that a modest change in the pre-stress conditions, from strike-slip to reverse stress field, strongly influences the magnitude of any plastic criterion and the extent of off-fault deformation. Using a Coulomb criterion, we also compute the local preferential orientations for failure and the corresponding type of faulting induced. We have shown that a simulation within the so-called “forbidden regime” will predict some reverse faulting in the off-fault medium.
Then, because of the various feedbacks that exist between the dynamic rupture and the bulk, we show that the discrepancy is even more pronounced when inelastic deformation can occur in the medium as the rupture propagates (figure 12). Both in plastic and dynamic damage models, the resulting pattern of inelastic deformation is significantly different. We show that a modest change in pre-stress conditions from initial reverse to initial strike-slip stress field would underestimate the damage zone width by a factor of 3 to 6 (Figure 12b,c). We would like to underline that, while previous numerical studies have modelled inelastic deformation under different stresses regimes (Shi & Ben-Zion, 2006; Dunham et al., 2011b; Templeton & Rice, 2008; Okubo et al., 2019), as illustrated in Figure 3a, they did not investigate the effect of pre-stress by keeping S ratio, stress drop, and angle Ψ simultaneously constant. By keeping these rupture parameters constant among our models, we demonstrate the importance of initial stress field only. We also show that the effect of far field stresses can have a significant impact on the rupture dynamics (Figure 10). In a passive, elastic medium, pre-stress indeed does not affect fault slip. In inelastic medium, such as displayed in Figures 9 and 11 the evolving medium through energy loss and/or trapped-waves influences back both the slip rate and rupture velocity on the fault. We observe that this effect is more predominant with higher amount of off-fault deformation, hence when the initial stress field is rightfully set up and/or if a pre-existing damage zone is modelled. Therefore in 2D plane-strain simulations, the initial three dimensional stress field is important to model an accurate evolution of the off-fault medium.

To conclude, pre-stresses can significantly affect both off-fault damage and on-fault rupture dynamics even if other key parameters are kept constant: cohesive zone, nucleation size, seismic ratio, stress drop. Although none of the presented numerical models are meant to reproduce field observations exactly the sheer increase in observations opens up the potential for statistical comparisons between models and observations. This makes it even more urgent to set-up the correct initial “3D” stress field even in “2D” numerical simulations.

## Acknowledgments

This study was supported by the Agence National de la Recherche (ANR) IDEAS contract ANR-19-CE31-0004-01. H. S. B. acknowledges the European Research Council grant PERSISMO (grant 865411) for partial support of this work. We thank Prof. Shiqing Xi and an anonymous reviewer for constructive feedback that greatly helped us in conveying our message. We would also like to thank our editor Prof. Fukuyama for his helpful insights in the process.

## Data Availability

This study uses numerical data only. All of the models and data sets are produced by the authors. The numerical model used to perform the simulations (Ampuero, 2012) is available at the following link https://github.com/jpampuero/sem2dpack. The module developed by Thomas et al. (2017b) to simulate dynamic damage will be shared on request to the corresponding author. We provide a Matlab code to check a given initial stress field and plot Figure 3 https://github.com/louisejeandetribes/PreStress_2Dplanestrain.git
References

Abdelmeguid, M., Ma, X., & Elbanna, A., 2019. A novel hybrid finite element-spectral boundary integral scheme for modeling earthquake cycles: Application to rate and state faults with low-velocity zones, *Journal of Geophysical Research: Solid Earth*, **124**(12), 12854–12881.

Ampuero, J., 2012. Sem2dpack, a spectral element software for 2d seismic wave propagation and earthquake source dynamics, v2. 3.8.

Anderson, E., 1951. The dynamics of faulting and dyke formation with applications to britain: Edinburgh, *Oliver and Boyd*, **206**.

Anderson, E. M., 1905. The dynamics of faulting, *Transactions of the Edinburgh Geological Society*, **8**(3), 387–402.

Andrews, D., 2005. Rupture dynamics with energy loss outside the slip zone, *Journal of Geophysical Research: Solid Earth*, **110**(B1).

Andrews, D. J., 1976. Rupture velocity of plane strain shear cracks, *J. Geophys. Res.*, **81**(B32), 5679–5689.

Barbot, S., Lapusta, N., & Avouac, J.-P., 2012. Under the hood of the earthquake machine: Toward predictive modeling of the seismic cycle, *Science*, **336**(6082), 707–710.

Ben-Zion, Y. & Rice, J. R., 1997. Dynamic simulations of slip on a smooth fault in an elastic solid, *Journal of Geophysical Research: Solid Earth*, **102**(B8), 17771–17784.

Ben-Zion, Y. & Shi, Z., 2005. Dynamic rupture on a material interface with spontaneous generation of plastic strain in the bulk, *Earth and Planetary Science Letters*, **236**(1-2), 486–496.

Ben-Zion, Y., Peng, Z., Okaya, D., Seeber, L., Armbruster, J. G., Ozer, N., Michael, A. J., Baris, S., & Aktar, M., 2003. A shallow fault-zone structure illuminated by trapped waves in the karadere-duzce branch of the north anatolian fault, western turkey, *Geophysical Journal International*, **152**(3), 699–717.

Bhat, H. S., Dmowska, R., Rice, J. R., & Kame, N., 2004. Dynamic slip transfer from the denali to totschunda faults, alaska: Testing theory for fault branching, *Bull. Seism. Soc. Am.*, **94**, S202–S213.

Bhat, H. S., Biegel, R. L., Rosakis, A. J., & Sammis, C. G., 2010. The effect of asymmetric damage on dynamic shear rupture propagation ii: With mismatch in bulk elasticity, *Tectonophysics*, **493**(3), 263–271.

Bhat, H. S., Rosakis, A. J., & Sammis, C. G., 2012a. A micromechanics based constitutive model for brittle failure at high strain rates, *J. Appl. Mech.*, **79**(3).

Bhat, H. S., Rosakis, A. J., & Sammis, C. G., 2012b. A micromechanics based constitutive model for brittle failure at high strain rates, *Journal of Applied Mechanics*, **79**(3).

Biegel, R. L. & Sammis, C. G., 2004. *Relating Fault Mechanics to Fault Zone Structure*, vol. 47, pp. 65–111, Elsevier.
Biegel, R. L., Bhat, H. S., Sammis, C. G., & Rosakis, A. J., 2010. The effect of asymmetric damage on dynamic shear rupture propagation i: No mismatch in bulk elasticity, *Tectonophysics*, **493**(3), 254–262.

Brenguier, F., Campillo, M., Hadziioannou, C., Shapiro, N. M., Nadeau, R. M., & Larose, E., 2008. Postseismic relaxation along the san andreas fault at parkfield from continuous seismological observations, *Science*, **321**(5895), 1478–1481.

Byerlee, J., 1978. Friction of rocks, in *Rock friction and earthquake prediction*, pp. 615–626, Springer.

Canitano, A., Godano, M., & Thomas, M. Y., 2021. Inherited state of stress as a key factor controlling slip and slip mode: inference from the study of a slow slip event in the longitudinal valley, taiwan, *Geophysical Research Letters*, **48**(3), e2020GL090278.

Cappa, F., Perrin, C., Manighetti, I., & Delor, E., 2014. Off-fault long-term damage: A condition to account for generic, triangular earthquake slip profiles, *Geochemistry, Geophysics, Geosystems*, **15**(4), 1476–1493.

Célérier, B., 2008. Seeking anderson's faulting in seismicity: a centennial celebration, *Reviews of Geophysics*, **46**(4).

Chester, F. M., Evans, J. P., & Biegel, R. L., 1993. Internal structure and weakening mechanisms of the san andreas fault, *Journal of Geophysical Research-Solid Earth*, **98**(B1), 771–786.

Cochran, E. S., Li, Y.-G., Shearer, P. M., Barbot, S., Fialko, Y., & Vidale, J. E., 2009. Seismic and geodetic evidence for extensive, long-lived fault damage zones, *Geology*, **37**(4), 315–318.

Collettini, C. & Sibson, R. H., 2001. Normal faults, normal friction?, *Geology*, **29**(10), 927–930.

Cowie, P. A. & Scholz, C. H., 1992. Growth of faults by accumulation of seismic slip, *Journal of Geophysical Research: Solid Earth*, **97**(B7), 11085–11095.

Cubas, N., Lapusta, N., Avouac, J.-P., & Perfettini, H., 2015. Numerical modeling of long-term earthquake sequences on the ne japan megathrust: Comparison with observations and implications for fault friction, *Earth and Planetary Science Letters*, **419**, 187–198.

Dalguer, L., Irikura, K., & Riera, J., 2003. Simulation of tensile crack generation by three-dimensional dynamic shear rupture propagation during an earthquake, *Journal of Geophysical Research: Solid Earth*, **108**(B3).

Das, S. & Aki, K., 1977. A numerical study of two-dimensional spontaneous rupture propagation, *Geophys. J. R. astr. Soc.*, **50**, 643–668.

Das, S. & Kostrov, B., 1986. Fracture of a single asperity on a finite fault: a model for weak earthquakes?, in *Earthquake Source Mechanics*, vol. 6 of *Geophys. Monograph*, pp. 91–96.

Das, S. & Kostrov, B. V., 1987. On the numerical boundary integral equation method for three-dimensional dynamic shear crack problems, *J. Appl. Mech.*, **54**(1), 99–104.
Day, S. M., Dalguer, L. A., Lapusta, N., & Liu, Y., 2005. Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture, *Journal of Geophysical Research-solid Earth*, **110**(B12), B12307.

Dor, O., Rockwell, T. K., & Ben-Zion, Y., 2006. Geological observations of damage asymmetry in the structure of the san jacinto, san andreas and punchbowl faults in southern california: A possible indicator for preferred rupture propagation direction, *pure and applied geophysics*, **163**(2), 301–349.

Duan, B. & Day, S. M., 2008. Inelastic strain distribution and seismic radiation from rupture of a fault kink, *Journal of Geophysical Research: Solid Earth*, **113**(B12).

Dunham, E. M., Belanger, D., Cong, L., & Kozdon, J. E., 2011a. Earthquake ruptures with strongly rate-weakening friction and off-fault plasticity, part 2: Nonplanar faults, *Bulletin of the Seismological Society of America*, **101**(5), 2308–2322.

Dunham, E. M., Belanger, D., Cong, L., & Kozdon, J. E., 2011b. Earthquake ruptures with strongly rate-weakening friction and off-fault plasticity, part 1: Planar faults, *Bulletin of the Seismological Society of America*, **101**(5), 2296–2307.

Erickson, B. A., Dunham, E. M., & Khosravifar, A., 2017. A finite difference method for off-fault plasticity throughout the earthquake cycle, *Journal of the Mechanics and Physics of Solids*, **109**, 50–77.

Erickson, B. A., Jiang, J., Barall, M., Lapusta, N., Dunham, E. M., Harris, R., Abrahams, L. S., Allison, K. L., Ampuero, J., Barbot, S., Cattania, C., Elbanna, A., Fialko, Y., Idini, B., Kozdon, J. E., Lambert, V., Liu, Y., Luo, Y., Ma, X., Best McKay, M., Segall, P., Shi, P., van den Ende, M., & Wei, M., 2020. The community code verification exercise for simulating sequences of earthquakes and aseismic slip (seas), *Seismological Research Letters*, **91**(2A), 874–890.

Faulkner, D., Mitchell, T., Healy, D., & Heap, M., 2006. Slip on'weak'faults by the rotation of regional stress in the fracture damage zone, *Nature*, **444**(7121), 922–925.

Faulkner, D. R., Mitchell, T. M., Jensen, E., & Cembrano, J., 2011. Scaling of fault damage zones with displacement and the implications for fault growth processes, *J. Geophys. Res.*, **116**(B5).

Fliss, S., Bhat, H. S., Dmowska, R., & Rice, J. R., 2005. Fault branching and rupture directivity, *J. Geophys. Res.*, **B06312**.

Froment, B., McGuire, J., van der Hilst, R., Gouédard, P., Roland, E., Zhang, H., & Collins, J., 2014. Imaging along-strike variations in mechanical properties of the gofar transform fault, east pacific rise, *Journal of Geophysical Research: Solid Earth*, **119**(9), 7175–7194.

Fukuyama, E., Ellsworth, W. L., Waldhauser, F., & Kubo, A., 2003. Detailed fault structure of the 2000 western tottori, japan, earthquake sequence, *Bulletin of the Seismological Society of America*, **93**(4), 1468–1478.

Gabriel, A. A., Ampuero, J.-P., Dalguer, L. A., & Mai, P. M., 2013. Source properties of dynamic rupture pulses with off-fault plasticity, *Journal of Geophysical Research-solid Earth*, **118**(8), 4117–4126.
Healy, D., Blenkinsop, T. G., Timms, N. E., Meredith, P. G., Mitchell, T. M., & Cooke, M. L., 2015. Polymodal faulting: Time for a new angle on shear failure, *Journal of Structural Geology*, **80**, 57–71.

Hiramatsu, Y., Honma, H., Saiga, A., Furumoto, M., & Ooida, T., 2005. Seismological evidence on characteristic time of crack healing in the shallow crust, *Geophysical Research Letters*, **32**(9).

Hok, S., Campillo, M., Cotton, F., Favreau, P., & Ionescu, I., 2010. Off-fault plasticity favors the arrest of dynamic ruptures on strength heterogeneity: Two-dimensional cases, *Geophysical Research Letters*, **37**, L02306.

Huang, Y., Ampuero, J.-P., & Helmberger, D. V., 2014. Earthquake ruptures modulated by waves in damaged fault zones, *J. Geophys. Res. Solid Earth*, **119**(4), 3133–3154.

Idini, B. & Ampuero, J. P., 2020. Fault-zone damage promotes pulse-like rupture and back-propagating fronts via quasi-static effects, *Geophysical Research Letters*, **47**(23), e2020GL090736.

Im, K., Saffer, D., Marone, C., & Avouac, J.-P., 2020. Slip-rate-dependent friction as a universal mechanism for slow slip events, *Nature Geoscience*, **13**(10), 705–710.

Jaeger, J., 1979. Cook., ngw fundamentals of rock mechanics, *Loydon, methuen*.

Johnson, S. E., Song, W. J., Vel, S. S., Song, B. R., & Gerbi, C. C., 2021. Energy partitioning, dynamic fragmentation, and off-fault damage in the earthquake source volume, *Journal of Geophysical Research: Solid Earth*, **126**(11), e2021JB022616, e2021JB022616 2021JB022616.

Johri, M., Dunham, E. M., Zoback, M. D., & Fang, Z., 2014. Predicting fault damage zones by modeling dynamic rupture propagation and comparison with field observations, *Journal of Geophysical Research: Solid Earth*, **119**(2), 1251–1272.

Kame, N., Rice, J. R., & Dmowska, R., 2003a. Effects of prestress state and rupture velocity on dynamic fault branching, *J. Geophys. Res.*, **108**(B5).

Kame, N., Rice, J. R., & Dmowska, R., 2003b. Effects of prestress state and rupture velocity on dynamic fault branching, *J. Geophys. Res.*, **108**(B5).

Kanamori, H., 2006. Lessons from the 2004 sumatra-andaman earthquake, *Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*, **364**(1845), 1927–1945.

Kaneko, Y., Ampuero, J. P., & Lapusta, N., 2011. Spectral-element simulations of long-term fault slip: Effect of low-rigidity layers on earthquake-cycle dynamics, *Journal of Geophysical Research: Solid Earth*, **116**(B10313).

Kato, N., 2004. Interaction of slip on asperities: Numerical simulation of seismic cycles on a two-dimensional planar fault with nonuniform frictional property, *Journal of Geophysical Research-Solid Earth*, **109**(B12), B12306.

Kelly, P., Sanderson, D., & Peacock, D., 1998. Linkage and evolution of conjugate strike-slip fault zones in limestones of somerset and northumbria, *Journal of Structural Geology*, **20**(11), 1477–1493.
King, G., Stein, R., & Lin, J., 1994. Coulomb static stress technique and the triggering of earthquakes, *Geophys. J. Int*, **146**, 747–759.

Kostrov, B. V., 1964. Selfsimilar problems of propagation of shear cracks, *J. Appl. Math. Mech.-USS.*, **28**(5), 1077–1087.

Lapusta, N., Rice, J. R., Ben-Zion, Y., & Zheng, G., 2000. Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate-and state-dependent friction, *Journal of Geophysical Research: Solid Earth*, **105**(B10), 23765–23789.

Lecomte, E., Jolivet, L., Lacombe, O., Denèle, Y., Labrousse, L., & Le Pourhiet, L., 2010. Geometry and kinematics of mykonos detachment, cyclades, greece: Evidence for slip at shallow dip, *Tectonics*, **29**(5), n/a–n/a.

Lisle, R. J., Orife, T. O., Arlegui, L., Liesa, C., & Srivastava, D. C., 2006. Favoured states of palaeostress in the earth's crust: evidence from fault-slip data, *Journal of Structural Geology*, **28**(6), 1051–1066.

Liu, D., Duan, B., & Luo, B., 2020. Eqsimu: a 3-d finite element dynamic earthquake simulator for multicycle dynamics of geometrically complex faults governed by rate- and state-dependent friction, *Geophysical Journal International*, **220**(1), 598–609.

Lyakhovsky, V., Ben-Zion, Y., & Agnon, A., 1997. Distributed damage, faulting, and friction, *Journal of Geophysical Research-solid Earth*, **102**(B12), 27635–27649.

Ma, S., 2008. A physical model for widespread near-surface and fault zone damage induced by earthquakes, *Geochem. Geophys. Geosyst.*, **9**(11).

Madariaga, R., 1976. Dynamics of an expanding circular fault, *Bull. Seism. Soc. Am.*, **66**(3), 639–666.

Manighetti, I., King, G., Gaudemer, Y., Scholz, C., & Doubre, C., 2001. Slip accumulation and lateral propagation of active normal faults in afar, *Journal of Geophysical Research: Solid Earth*, **106**(B7), 13667–13696.

Manighetti, I., King, G., & Sammis, C. G., 2004. The role of off-fault damage in the evolution of normal faults, *Earth and Planetary Science Letters*, **217**(3-4), 399–408.

Mitchell, T. & Faulkner, D., 2009. The nature and origin of off-fault damage surrounding strike-slip fault zones with a wide range of displacements: A field study from the atacama fault system, northern chile, *Journal of Structural Geology*, **31**(8), 802–816.

Ngo, D., Huang, Y., Rosakis, A., Griffith, W., & Pollard, D., 2012. Off-fault tensile cracks: A link between geological fault observations, lab experiments, and dynamic rupture models, *Journal of Geophysical Research: Solid Earth*, **117**(B1).

Okubo, K., Bhat, H. S., Rougier, E., Marty, S., Schubnel, A., Lei, Z., Knight, E. E., & Klinger, Y., 2019. Dynamics, radiation, and overall energy budget of earthquake rupture with coseismic off-fault damage, *Journal of Geophysical Research: Solid Earth*, **124**(11), 11771–11801.
Ozawa, S. & Ando, R., 2021. Mainshock and aftershock sequence simulation in geometrically complex fault zones, *J. Geophys. Res.*, **126**(2).

Palmer, A. C. & Rice, J. R., 1973. The growth of slip surfaces in the progressive failure of over-consolidated clay, *Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences*, **332**(1591), 527-548.

Poliakov, A. N., Dmowska, R., & Rice, J. R., 2002. Dynamic shear rupture interactions with fault bends and off-axis secondary faulting, *Journal of Geophysical Research: Solid Earth*, **107**(B11), ESE-6.

Preuss, S., Ampuero, J. P., Gerya, T. V., & van Dinther, Y., 2020. Characteristics of earthquake ruptures and dynamic off-fault deformation on propagating faults, *Solid Earth*.

Rice, J. R., 2002. New perspectives on crack and fault dynamics, in *Mechanics for a New Millennium: Proceedings of the 20th International Congress of Theoretical and Applied Mechanics Chicago, Illinois, USA 27 August - 2 September 2000*, pp. 1-24, Springer Netherlands, Dordrecht.

Rice, J. R., Sammis, C. G., & Parsons, R., 2005. Off-fault secondary failure induced by a dynamic slip pulse, *Bulletin of the Seismological Society of America*, **95**(1), 109-134.

Richards-Dinger, K. B. & Shearer, P. M., 2000. Earthquake locations in southern California obtained using source-specific station terms, *Journal of Geophysical Research-Solid Earth*, **105**(B5), 10939-10960.

Rodriguez Padilla, A. M., Oskin, M. E., Milliner, C. W. D., & Plesch, A., 2022. Accrual of widespread rock damage from the 2019 ridgecrest earthquakes, *Nature Geoscience*, **15**(3), 222-226.

Romanet, P., Bhat, H. S., Jolivet, R., & Madariaga, R., 2018. Fast and slow slip events emerge due to fault geometrical complexity, *Geophysical Research Letters*, **45**(10), 4809-4819.

Sammis, C. G., Rosakis, A. J., & Bhat, H. S., 2009. Effects of off-fault damage on earthquake rupture propagation: Experimental studies, *Pure Appl. Geophys.*, **166**.

Savage, H. M. & Brodsky, E. E., 2011. Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands in fault damage zones.

Shi, Z. & Ben-Zion, Y., 2006. Dynamic rupture on a bimaterial interface governed by slip-weakening friction, *Geophysical Journal International*, **165**(2), 469-484.

Shipton, Z. & Cowie, P., 2001. Damage zone and slip-surface evolution over µm to km scales in high-porosity navajo sandstone, utah, *Journal of Structural Geology*, **23**(12), 1825-1844.

Sibson, R., Ghisetti, F., & Crookbain, R., 2012. Andersonian wrench faulting in a regional stress field during the 2010-2011 canterbury, new zealand, earthquake sequence, *Geological Society, London, Special Publications*, **367**(1), 7-18.

Stein, R. S., Barka, A. A., & Dieterich, J. H., 1997. Progressive failure on the north anatolian fault since 1939 by earthquake stress triggering, *Geophysical Journal International*, **128**(3), 594-604.
Tal, Y. & Faulkner, D. R., 2022. The effect of fault roughness and earthquake ruptures on the evolution and scaling of fault damage zones, *Journal of Geophysical Research: Solid Earth*, **127**.

Templeton, E. L. & Rice, J. R., 2008. Off-fault plasticity and earthquake rupture dynamics: 1. dry materials or neglect of fluid pressure changes, *Journal of Geophysical Research: Solid Earth*, **113**(B9).

Thakur, P. & Huang, Y., 2021. Influence of fault zone maturity on fully dynamic earthquake cycles, *Geophysical Research Letters*, **48**(17), e2021GL094679.

Thomas, M. Y. & Bhat, H. S., 2018. Dynamic evolution of off-fault medium during an earthquake: a micromechanics based model, *Geophysical Journal International*, **214**(2), 1267–1280.

Thomas, M. Y., Lapusta, N., Noda, H., & Avouac, J.-P., 2014. Quasi-dynamic versus fully dynamic simulations of earthquakes and aseismic slip with and without enhanced coseismic weakening, *Journal of Geophysical Research: Solid Earth*, **119**(3), 1986–2004.

Thomas, M. Y., Avouac, J.-P., & Lapusta, N., 2017a. Rate-and-state friction properties of the longitudinal valley fault from kinematic and dynamic modeling of seismic and aseismic slip, *Journal of Geophysical Research: Solid Earth*, **122**(4), 3115–3137.

Thomas, M. Y., Bhat, H. S., & Klinger, Y., 2017b. Effect of brittle off-fault damage on earthquake rupture dynamics, *Fault zone dynamic processes: Evolution of fault properties during seismic rupture*, **227**, 255.

Townend, J. & Zoback, M. D., 2000. How faulting keeps the crust strong, *Geology*, **28**(5), 399–402.

Tullis, T. E. & Schubert, G., 2015. 4.06 - mechanisms for friction of rock at earthquake slip rates, in *Treatise on Geophysics (Second Edition)*, pp. 139–159, Elsevier, Oxford.

Uphoff, C., May, D. A., & Gabriel, A.-A., 2022. A discontinuous Galerkin method for sequences of earthquakes and aseismic slip on multiple faults using unstructured curvilinear grids, *Geophysical Journal International*, ggac467.

van den Ende, M. P. A., Chen, J., Ampuero, J. P., & Niemeijer, A. R., 2018. A comparison between rate-and-state friction and microphysical models, based on numerical simulations of fault slip, *Tectonophysics*, **733**, 273–295.

Vermilye, J. M. & Scholz, C. H., 1998. The process zone: A microstructural view of fault growth, *Journal of Geophysical Research-solid Earth*, **103**(B6), 12223–12237.

Walsh, J. B., 1965a. The effect of cracks in rocks on poisson’s ratio, *J. Geophys. Res.*, **70**(20), 5249–5257.

Walsh, J. B., 1965b. The effect of cracks on the compressibility of rock, *J. Geophys. Res.*, **70**(2), 381–389.

Wilson, J. E., Chester, J. S., & Chester, F. M., 2003. Microfracture analysis of fault growth and wear processes, punchbowl fault, san andreas system, california, *Journal of Structural Geology*, **25**(11), 1855–1873.
Xu, S., Ben-Zion, Y., Ampuero, J.-P., & Lyakhovsky, V., 2014. Dynamic ruptures on a frictional interface with off-fault brittle damage: Feedback mechanisms and effects on slip and near-fault motion, *Pure and Applied Geophysics*, **172**(5), 1243–1267.

Yamashita, T., 2000. Generation of microcracks by dynamic shear rupture and its effects on rupture growth and elastic wave radiation, *Geophysical Journal International*, **143**(2), 395–406.

Yukutake, Y., Iio, Y., Katao, H., & Shibutani, T., 2007. Estimation of the stress field in the region of the 2000 western tottori earthquake: Using numerous aftershock focal mechanisms, *Journal of Geophysical Research: Solid Earth*, **112**(B9).

Zoback, M. D., Zoback, M. L., Mount, V. S., Suppe, J., Eaton, J. P., Healy, J. H., Oppenheimer, D., Reasenberg, P., Jones, L., Raleigh, C. B., Wong, I. G., Scotti, O., & Wentworth, C., 1987. New evidence on the state of stress of the san andreas fault system, *Science*, **238**(4830), 1105–1111.
Supplementary Figures

S1 Initial stress field

This section provides the Mohr-Coulomb diagrams that illustrate the initial stress fields of the four
main models of the manuscript.

S2 Elastic models - dynamic stress changes

This section displays the Drucker-Prager criterion for the model with $\Psi = 55^\circ$ and $S = 2$, computed
with the whole stress tensor and the in-plane stresses only.

S3 Supplementary models with plastic deformation

This section provides the supplementary figure cited in section 4.2.1: *Role of the initial stress field on plastic deformation.*
<figure><img src="image_21.png" /><figcaption>Figure S1: Initial stress fields of the four main models. The Mohr Circle is depicted between $σ_1$ and $σ_2$ when the initial stress field is incorrect (i.e., it actually corresponds to a reverse stress field in 3D), and between $σ_1$ and $σ_3$ when the initial stress field is correct (i.e., it does correspond to a strike-slip stress field in 3D). It is important to notice that we artificially trigger the rupture by applying a higher shear stress on a nucleation patch that is not included in this figure.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure S2: Drucker-Prager criterion for the model with $\Psi = 55^\circ$ and $S = 2$, computed with the whole stress tensor (a) and the in-plane stresses only (b). Contrary to Figure 5, here we only show the dynamic change of plastic criterion by removing the contribution from the initial stress field ($F_{DP} - F_{DP,ini}$). Be careful, here the simulations are identical. Only the error in computing the Drucker-Prager criteria, by solely using the in-plane stresses for case (b), is responsible for the discrepancy.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure S3: (a) & (b), plastic deformation for the models with $\Psi = 55^{\circ}$, S=2, with initial strike-slip stress field (case i) and reverse stress field (case iv). (c) & (d), plastic deformation for the models with $\Psi = 15^{\circ}$, S=1, with initial strike-slip stress field (case vi) and reverse stress field (case xii). (e) & (f), plastic deformation for the models with $\Psi = 30^{\circ}$, S=1, with initial strike-slip stress field (case xiii) and reverse stress field (casexiv). In all cases the cohesion $c = 0$. (g) Criterion for accurate initial stress field under plane-strain approximation displayed as $\gamma - \mu_0$. Red circles give the initial parameters for the four simulations.</figcaption></figure>