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

Manuscript submitted to Geophys. J. Int.

Louise Jeandet Ribes1, Marion Y. Thomas1 and Harsha S. Bhat2

  1. 1.

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

  2. 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
mediums, 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 the 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, friction, numerical modelling

##  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 the 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, 1976b; 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 year (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 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 order, 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, 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, people 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 have 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 an
homogenized damage mechanics based constitutive laws (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 others.

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 constrains 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 off-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
seismicaly 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 tractions acts.
Those planes are perpendicular one to another and are called the principal
planes. The corresponding stresses acting along the three principal directions
are called the principal stresses with the following inequality (stresses
positive in tension)

$\sigma_{1}<\sigma_{2}<\sigma_{3}$ (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}$, this Coulomb failure criterion may be written as:

$\tau_{y}=c+\mu\sigma_{eff}=c+\mu(\sigma_{n}-P_{f})$ (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$ correspond to the slope of the failure envelope on
a Mohr diagram and $\theta$ corresponds to the angle between the normal to the
fault plane and $\sigma_{1}$ (Figure 1). Hence we can write:

$\begin{array}[]{ll}\sigma_{yy}=&\left(\dfrac{\sigma_{1}+\sigma_{3}}{2}\right)-\left(\dfrac{\sigma_{3}-\sigma_{1}}{2}\right)\cos{2\theta}\\\
\sigma_{xx}=&\left(\dfrac{\sigma_{1}+\sigma_{3}}{2}\right)+\left(\dfrac{\sigma_{3}-\sigma_{1}}{2}\right)\cos{2\theta}\\\
\sigma_{xy}=&\left(\dfrac{\sigma_{3}-\sigma_{1}}{2}\right)\sin{2\theta}\\\
\end{array}$ (3)

Figure 1: a) Mohr–Coulomb criterion in ($\sigma,\tau$)-space. b), c) & d)
Orientation of the failure plane relative to the largest principal stress for
normal, strike-slip, and reverse faulting respectively.

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

$R=\left(\dfrac{\sigma_{3}-\sigma_{1}}{2}\right)=\sqrt{\left(\dfrac{\sigma_{xx}-\sigma_{yy}}{2}\right)^{2}+\sigma_{xy}^{2}}$
(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{array}[]{ccc}|\sigma_{yy}|<|\sigma_{xx}|&\mbox{if}&0\degree<\Psi\leq
45\degree\\\ |\sigma_{yy}|\geq|\sigma_{xx}|&\mbox{if}&45\degree<\Psi\leq
90\degree\end{array}$ (5)

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 envelop 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$\degree$ and 68$\degree$ and
the new fault should make an angle $\Psi$ with $\sigma_{1}$ between
32$\degree$ and 22$\degree$ 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), fault are optimally oriented when they makes an angle
$\Psi\simeq 30\degree$ with $\sigma_{1}$.

Anderson theory is based on the assumption that Earth’s free surface requires
a principal stress direction to be subvertical ($\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\degree$ and
$\omega\approx 60\degree$ 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$\degree$ 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
Darfield earthquake Mw 7.1 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. Hence, including
these exceptions, it suggests that the Anderson 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 in prescribed in
numerical models

In 2D 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 $\Psi$
(the angle between $\sigma_{1}$ and the primary fault). Steep $\Psi$ favors
inelasticity on the extensional side and shallow $\Psi$ on the compressional
side. Moreover initial stresses seems 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 $\sigma_{1}$ and $\sigma_{3}$ and an angle
$\Psi$ 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 $\sigma_{1}$ and $\sigma_{3}$, then the
fault is in a proper strike-slip stress field. If the two in-plane stresses
are $\sigma_{1}$ and $\sigma_{2}$, 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 greater 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, $\sigma_{ij}$, (tensile positive) is
represented as,

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

For convenience, we define the following ratios:

$\gamma=\dfrac{\sigma_{xx}}{\sigma_{yy}}\textrm{ ~{}\&~{}
}\mu_{0}=\dfrac{\sigma_{xy}}{-\sigma_{yy}}$ (7)

which leads to:

$\tan{2\Psi}=-\tan{2\theta}=-\dfrac{2\sigma_{xy}}{\sigma_{xx}-\sigma_{yy}}=\dfrac{2\mu_{0}}{\gamma-1}\\\
$ (8)

Let $\sigma_{1},\sigma_{2}$ and $\sigma_{3}$ represent the maximum,
intermediate and minimum compressive principal stresses. For a strike-slip
fault, $\sigma_{1}$ and $\sigma_{3}$ should lie on the $x-y$ plane and
$\sigma_{2}$ should be parallel to the $z$ axis, i.e out of plane for a 2D
simulation (Figure 1). Hence,

$\begin{array}[]{l}\sigma_{1}=\left(\dfrac{\sigma_{xx}+\sigma_{yy}}{2}\right)-R=\left(\dfrac{\sigma_{xx}+\sigma_{yy}}{2}\right)\left(1+\sqrt{\left(\dfrac{\gamma-1}{\gamma+1}\right)^{2}+\dfrac{4\mu_{0}^{2}}{(\gamma+1)^{2}}}\right),\\\
\sigma_{3}=\left(\dfrac{\sigma_{xx}+\sigma_{yy}}{2}\right)+R=\left(\dfrac{\sigma_{xx}+\sigma_{yy}}{2}\right)\left(1-\sqrt{\left(\dfrac{\gamma-1}{\gamma+1}\right)^{2}+\dfrac{4\mu_{0}^{2}}{(\gamma+1)^{2}}}\right).\\\
\sigma_{2}=\sigma_{zz}=\rho gz(1-\lambda)\\\ \end{array}$ (9)

where $\lambda$ is the pore pressure coefficient. 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)$
(10)

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

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}\left(-\sigma_{yy}\right)$ over a length $L_{c}$.

We require that both $\sigma_{1}$ and $\sigma_{3}$ should be compressive, i.
e. $\sigma_{1}<\sigma_{3}<0$. This implies that

${\sqrt{\left(\dfrac{\gamma-1}{\gamma+1}\right)^{2}+\dfrac{4\mu_{0}^{2}}{(\gamma+1)^{2}}}~{}<~{}1}$
(11)

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

$\begin{array}[]{ll}\sigma_{1}<\sigma_{2}&\Rightarrow\left(\dfrac{\sigma_{xx}+\sigma_{yy}}{2}\right)\left(1+\sqrt{\left(\dfrac{\gamma-1}{\gamma+1}\right)^{2}+\dfrac{4\mu_{0}^{2}}{(\gamma+1)^{2}}}\right)~{}<~{}\nu(\sigma_{xx}+\sigma_{yy})\\\
&\Rightarrow{\sqrt{\left(\dfrac{\gamma-1}{\gamma+1}\right)^{2}+\dfrac{4\mu_{0}^{2}}{(\gamma+1)^{2}}}~{}>~{}2\nu-1}\end{array}$
(12)
$\begin{array}[]{ll}\sigma_{3}>\sigma_{2}&\Rightarrow\left(\dfrac{\sigma_{xx}+\sigma_{yy}}{2}\right)\left(1-\sqrt{\left(\dfrac{\gamma-1}{\gamma+1}\right)^{2}+\dfrac{4\mu_{0}^{2}}{(\gamma+1)^{2}}}\right)~{}>~{}\nu(\sigma_{xx}+\sigma_{yy})\\\
&\Rightarrow{\sqrt{\left(\dfrac{\gamma-1}{\gamma+1}\right)^{2}+\dfrac{4\mu_{0}^{2}}{(\gamma+1)^{2}}}~{}>~{}1-2\nu}\end{array}$
(13)

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

${1-2\nu~{}<~{}\sqrt{\left(\dfrac{\gamma-1}{\gamma+1}\right)^{2}+\dfrac{4\mu_{0}^{2}}{(\gamma+1)^{2}}}~{}<~{}1}$
(14)

This inequality is plotted in Figure 3 as a function of $\gamma$ and
$\mu_{0}$. The area in grey defines 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 3b, we plot the same criterion as a function of $\Psi$ and the
seismic ratio $S$, as defined by (Andrews, 1976a; Das & Aki, 1977),

$S=\dfrac{\mu_{s}(-\sigma_{yy})-\sigma_{xy}}{\sigma_{xy}-\mu_{d}(-\sigma_{yy})}=\dfrac{\mu_{s}-\mu_{0}}{\mu_{0}-\mu_{d}}$
(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 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.

Figure 3: Criterion for accurate initial stress field under plane-strain
approximation displayed as $\gamma-\mu_{0}$ plot (a) and $\Psi-S$ plot (b).
The gray areas represents the conditions for which the initial stresses favors
a strike-slip motion, i.e. when equation (14) is satisfied. In the
$\gamma-\mu_{0}$ plot, this criterion 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. Shaded red areas shows the two sets of parameters explored in
this study, with $\Psi=15^{o}$ and $\Psi=55^{o}$. The red empty circles
correspond to the end members. In subplot (b) 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).

##  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 pre-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
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}}$ (16)

where $G$ is the fracture energy, defined as:

$G=\frac{1}{2}D_{c}(\tau_{p}-\tau_{r})$ (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})}$
(18)

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

###  3.2 Initial parameters

In order to compare our study to the available literature, we run two sets of
models with $\Psi=55\degree$ and $\Psi=15\degree$ 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\degree$ and cases vii, viii, ix for
$\Psi=15\degree$) and three within the so-called forbidden regime (cases iv,
v, vi for $\Psi=55\degree$ and cases x, xi, xii for $\Psi=15\degree$). 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).

Table 1: Constants used in all models Parameter | Symbol | Value
---|---|---
depth | $z$ | -2.5 km
material density | $\rho$ | 2700 kg.m-3
S-wave velocity | $c_{s}$ | 3120 km.s-1
P-wave velocity | $c_{p}$ | 5600 km.s-1
Poisson’s ratio | $\nu$ | 0.27
characteristic slip | $D_{c}$ | 1m
pore pressure coefficient | $\lambda$ | 0.4

Then, to set up the remaining parameters we adopt the following strategy.
First we set the S ratio (equation 15), favouring a value that would lead to
subshear rupture ($S=2$) for $\Psi=55\degree$. However, for $\Psi=15$, 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$ (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})$ (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 gz(1-\lambda)$ (20)

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

$\begin{array}[]{ll}\text{eq.~{}\ref{eqn:ratio} \&
\ref{eqn:out_of_plane_stress}}&\Rightarrow\gamma=\dfrac{\sigma_{zz}}{\nu\sigma_{yy}}-1\\\\[10.0pt]
\text{eq.~{}\ref{eqn:stressdrop} \&
\ref{eqn:Sratio}}&\Rightarrow\gamma=\dfrac{\sigma_{zz}\left(\mu_{0}-\mu_{s}\right)}{\nu
S\Delta\tau}-1\\\\[10.0pt]
\text{eq.~{}\ref{eqn:Psi}}&\Rightarrow\gamma=\dfrac{\sigma_{zz}\left(0.5\left(\gamma-1\right)\tan{2\Psi}-\mu_{s}\right)}{\nu
S\Delta\tau}-1\\\\[10.0pt]
&\Rightarrow\gamma=\dfrac{\left(\alpha+\beta\right)+\mu_{s}\sigma_{zz}}{\left(\beta-\alpha\right)}\\\\[10.0pt]
&\text{with }\alpha=\nu S\Delta\tau\text{ \&
}\beta=0.5\sigma_{zz}\tan{2\Psi}\end{array}$ (21)

Knowing $\gamma$, we can derive:

$\mu_{0}=\dfrac{\left(\gamma-1\right)}{2}\tan{2\Psi}$ (22)
$\mu_{d}=\mu_{0}-\dfrac{\mu_{s}-\mu_{0}}{S}$ (23)
$\sigma_{yy}=\dfrac{\sigma_{zz}}{\nu\left(\gamma-1\right)}$ (24)
$\sigma_{xx}=\gamma\sigma_{yy}$ (25) $\sigma_{xy}=\mu_{0}(-\sigma_{yy})$ (26)

The values of initials parameters for the end-members models are summarized in
Table 2. 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.

| Parameter | Symbol | Initial stress field
---|---|---|---
| case (i) | case (vi) | case (vii) | case (xii)
|  |  | strike-slip | reverse | strike-slip | reverse
Input | angle | $\Psi$ | $55^{o}$ | $55^{o}$ | $15^{o}$ | $15^{o}$
seismic ratio | $S$ | 2 | 2 | 1 | 1
stress drop (MPa) | $\Delta\sigma$ | 10.2 | 10.2 | 10.6 | 10.6
static friction | $\mu_{s}$ | 0.68 | 0.52 | 0.68 | 0.52
Resulting parameters | dynamic friction | $\mu_{d}$ | 0.34 | 0.16 | 0.19 | 0.09
full stress tensor (MPa) | $\sigma_{yy}$ | -91 | -84 | -43 | -50
| $\sigma_{xx}$ | -61 | -67 | -108 | -102
| $\sigma_{zz}$ | -42 | -42 | -42 | -42
| $\sigma_{xy}$ | 41 | 23 | 19 | 15
principal stresses (MPa) | $\sigma_{1}$ | -120 | -101 | -113 | -106
| $\sigma_{2}$ | -42 | -51 | -42 | -45
| $\sigma_{3}$ | -32 | -42 | -38 | -42
process zone (m) | $R_{0}$ | 1047 | 1047 | 1570 | 1570
| nucleation length (m) | $L_{c}$ | 6562 | 6562 | 4374 | 4374
Table 2: Initial parameters used for the two end-members of each set of
models.

##  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.

#### 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$ (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}$ (28)

To illustrate the contribution of the out-of-plane stress in computing the
invariants we use case (i) with $\Psi=55^{o}$ 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 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 $\sim$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 quadrant. 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 is higher. Hence this leads to an underestimation of
the plastic deformation. Changing the initial state of stress, using case
(vii), i.e. for $\Psi=15^{o}$ 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.

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 $\Psi$ = 55 $\degree$. 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 $\Delta\tau$ (equation 19). Slip rate on the fault
(black curves) is super-imposed on the graphs. Figure 5:  Drucker-Prager
criterion for case (i) with $S=2$ and $\Psi$ = 55 $\degree$ (a,b) and for case
(vii) with $S=1$ and $\Psi=15^{o}$ (d,e), at t = 2.2 seconds. We compute the
invariants accounting for $S_{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. 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.

#### 4.1.2 Role of the initial stress field for 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^{o}$, $S$ = 2 and for $\Psi=15^{o}$, $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 a 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}$.

For initial reverse stress conditions, when $\Psi=55^{o}$ (Figure 6b, case
vi), positive $F_{DP}$ is observed in the tensile quadrants. When
$\Psi=15^{o}$ (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. Figure 7 shows the continuum of all models for
$\Psi=55^{o}$. 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 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.

Figure 6: Drucker-Prager criterion for cases (i) & (vi), with $S=2$ and $\Psi$
= 55 $\degree$ (a,b) and for cases (vii) & (xii), with $S=1$ and $\Psi=15^{o}$
(d,e). Cases (i) & (vii) corresponds 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. Figure 7:
Drucker-Prager criterion for all models with $\Psi=55^{o}$ and S=2. We plot
the maximum value $F_{DP}$ between t=0 and t=4 seconds. Only positive values
are shown. 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 than 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.

#### 4.1.3 Role of the initial stress field for off fault rupture modes

We have shown that pre-stress conditions influence the values of $F_{DP}$
criterion induced by the rupture. Yet the Durcker-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\sigma_{xy}+\mu_{s}\Delta\sigma_{yy})$ (29)

where $\Delta\sigma_{xy}$ and $\Delta\sigma_{yy}$ are the change of shear and
normal stress respectively, 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 than for Figure 6, i.e., the two end-members for
$\Psi=55^{o}$, S = 2 (cases i & vi) and for $\Psi=15^{o}$, 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^{o}$. Likewise, the areas likely to
induced off-fault deformation are larger when the initial stress-field is
properly set up, and this for the two tested values of $\Psi$.

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 is different (Figure 8c). 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).

Figure 8: Coulomb stress change $\Delta CFF$ due to the dynamic rupture,
computed on optimally oriented planes and normalized by the dynamic stress
drop $\Delta\tau$ (equation 19), for cases (i) & (vi), with $S=2$ and $\Psi$ =
55 $\degree$ (a,b) and for cases (vii) & (xii), with $S=1$ and $\Psi=15^{o}$
(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 $\Delta 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 $\Delta CFF>0.2$. Slip rate on the
fault (grey curve) is super-imposed on the graphs.

###  4.2 Dynamic simulation with inelastic rheology

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

#### 4.2.1 Role of the initial stress field for 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. 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).

Figure 9: Plastic deformation for the models with $\Psi=55^{o}$, 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.

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

$\gamma=\sqrt{2e_{ij}e_{ij}}$ (30)

where $e_{ij}=e_{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 small 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 S1). For initial strike-slip stress conditions, the rupture
decelerates rapidly, its propagation prevented by the intense inelastic
deformation of the medium. The off-fault deformation is localised and
optimally oriented with respect to the far field stresses 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^{o}$ and $S$=1 Figure S2), little
deformation occurs for the strike-slip case (vii) and none for the reverse
case (xii). When $\Psi=30^{o}$ (optimally oriented fault, Figure S3),
deformation only occurs on the main fault. Therefore, a flat fault has to be
mis-oriented to produce significant off-fault deformation.

Figure 10: Cumulative slip (a-c-e-g) and slip rate (b-d-f-h) plotted every 0.4
seconds for $\Psi=55^{o}$ and $S=2$, with initial reverse stress fields (dark
solid lines, case vi) and initial strike-slip stress field (dotted color
lines, case i).

#### 4.2.2 Role of the initial stress field for dynamic damage

The last set of simulations use a micromechanics-based model to determine the
dynamically-triggered off-fault damage and its counter-impact on the rupture
dynamic. Inelastic deformation can occur in the model by either cracks opening
or cracks 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. The current inelastic state of the medium is
defined by the scalar $D$, the fraction of volume occupied by microcracks:

$D=\dfrac{4\pi}{3}N_{v}\left(a\cos\Phi_{c}+l\right)^{3}$ (31)

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$) 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. $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-members
cases (i) and (vi) with $\Psi=55^{o}$ 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. 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 above, 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 11: Simulation of a dynamic rupture with off-fault damage within an
homogenous 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 $\Psi$ = 55 $\degree$: 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.

##  5 Conclusion

Figure 12: Synthesis of the effect of initial stress on damage zone width for
cases (i) and (vi) with $\Psi=55^{o}$ and S=2, when the rupture has reached
seven times the 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 an
uniform initial damage density $D_{0}=0.1$). The shaded areas are computed by
changing by $\pm 10\%$ the threshold value.

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 principal stresses, at constant angle $\Psi$, have 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 small 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 feedbacks that exists 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 small 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 $\Psi$ 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.

###### Acknowledgements.

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.

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 here.

## References

  * Abdelmeguid et al. (2019) 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 (2012) Ampuero, J., 2012.  Sem2dpack, a spectral element software for 2d seismic wave propagation and earthquake source dynamics, v2. 3.8.
  * Anderson (1951) Anderson, E., 1951.  The dynamics of faulting and dyke formation with applications to britain: Edinburgh, Oliver and Boyd, 206.
  * Anderson (1905) Anderson, E. M., 1905.  The dynamics of faulting, Transactions of the Edinburgh Geological Society, 8(3), 387–402.
  * Andrews (1976a) Andrews, D., 1976a.  Rupture velocity of plane strain shear cracks, Journal of Geophysical Research, 81(32), 5679–5687.
  * Andrews (2005) Andrews, D., 2005.  Rupture dynamics with energy loss outside the slip zone, Journal of Geophysical Research: Solid Earth, 110(B1).
  * Andrews (1976b) Andrews, D. J., 1976b.  Rupture velocity of plane strain shear cracks, J. Geophys. Res., 81(B32), 5679–5689.
  * Barbot et al. (2012) 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 & Rice (1997) 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 & Shi (2005) 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 et al. (2003) 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 et al. (2004) 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 et al. (2010) 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 et al. (2012a) 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 et al. (2012b) 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 & Sammis (2004) Biegel, R. L. & Sammis, C. G., 2004.  Relating Fault Mechanics to Fault Zone Structure, vol. 47, pp. 65–111, Elsevier.
  * Biegel et al. (2010) 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 et al. (2008) 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 (1978) Byerlee, J., 1978.  Friction of rocks, in Rock friction and earthquake prediction, pp. 615–626, Springer.
  * Canitano et al. (2021) 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 et al. (2014) 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 (2008) Célérier, B., 2008.  Seeking anderson’s faulting in seismicity: a centennial celebration, Reviews of Geophysics, 46(4).
  * Chester et al. (1993) 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 et al. (2009) 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 & Sibson (2001) Collettini, C. & Sibson, R. H., 2001.  Normal faults, normal friction?, Geology, 29(10), 927–930.
  * Cowie & Scholz (1992) 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 et al. (2015) 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 et al. (2003) 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 & Aki (1977) Das, S. & Aki, K., 1977.  A numerical study of two-dimensional spontaneous rupture propagation, Geophys. J. R. astr. Soc., 50, 643–668.
  * Das & Kostrov (1986) 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 & Kostrov (1987) 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 et al. (2005) 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 et al. (2006) 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 & Day (2008) 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 et al. (2011a) 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 et al. (2011b) 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 et al. (2017) 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 et al. (2020) 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 et al. (2006) 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 et al. (2011) 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 et al. (2005) Fliss, S., Bhat, H. S., Dmowska, R., & Rice, J. R., 2005.  Fault branching and rupture directivity, J. Geophys. Res., B06312.
  * Froment et al. (2014) 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 et al. (2003) 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 et al. (2013) 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.
  * Hiramatsu et al. (2005) 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 et al. (2010) 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 et al. (2014) 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 & Ampuero (2020) 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 et al. (2020) 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 (1979) Jaeger, J., 1979.  Cook., ngw fundamentals of rock mechanics, Loydon, methuen.
  * Johnson et al. (2021) 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 et al. (2014) 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 et al. (2003a) 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 et al. (2003b) 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 (2006) 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 et al. (2011) 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 (2004) 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 et al. (1998) 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 et al. (1994) King, G., Stein, R., & Lin, J., 1994.  Coulomb static stress technique and the triggering of earthquakes, Geophys. J. Int, 146, 747–759.
  * Kostrov (1964) Kostrov, B. V., 1964.  Selfsimilar problems of propagation of shear cracks, J. Appl. Math. Mech.-USS., 28(5), 1077–1087.
  * Lapusta et al. (2000) 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 et al. (2010) 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 et al. (2006) 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 et al. (2020) 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 et al. (1997) Lyakhovsky, V., Ben-Zion, Y., & Agnon, A., 1997.  Distributed damage, faulting, and friction, Journal of Geophysical Research-solid Earth, 102(B12), 27635–27649.
  * Ma (2008) Ma, S., 2008.  A physical model for widespread near-surface and fault zone damage induced by earthquakes, Geochem. Geophys. Geosyst., 9(11).
  * Madariaga (1976) Madariaga, R., 1976.  Dynamics of an expanding circular fault, Bull. Seism. Soc. Am., 66(3), 639–666.
  * Manighetti et al. (2001) 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 et al. (2004) 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 & Faulkner (2009) 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 et al. (2012) 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 et al. (2019) 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 & Ando (2021) Ozawa, S. & Ando, R., 2021.  Mainshock and aftershock sequence simulation in geometrically complex fault zones, J. Geophys. Res., 126(2).
  * Palmer & Rice (1973) 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 et al. (2002) 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 et al. (2020) 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 (2002) 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 et al. (2005) 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 & Shearer (2000) 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 et al. (2022) 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 et al. (2018) 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 et al. (2009) 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 & Brodsky (2011) Savage, H. M. & Brodsky, E. E., 2011.  Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands in fault damage zones.
  * Shi & Ben-Zion (2006) 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 & Cowie (2001) Shipton, Z. & Cowie, P., 2001.  Damage zone and slip-surface evolution over $\mu$m to km scales in high-porosity navajo sandstone, utah, Journal of Structural Geology, 23(12), 1825–1844.
  * Sibson et al. (2012) 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 et al. (1997) 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 & Faulkner (2022) 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 & Rice (2008) 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 & Huang (2021) Thakur, P. & Huang, Y., 2021.  Influence of fault zone maturity on fully dynamic earthquake cycles, Geophysical Research Letters, 48(17), e2021GL094679.
  * Thomas & Bhat (2018) 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 et al. (2014) 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 et al. (2017a) 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 et al. (2017b) 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 & Zoback (2000) Townend, J. & Zoback, M. D., 2000.  How faulting keeps the crust strong, Geology, 28(5), 399–402.
  * Tullis & Schubert (2015) 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 et al. (2022) 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 et al. (2018) 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 & Scholz (1998) 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 (1965a) Walsh, J. B., 1965a.  The effect of cracks in rocks on poisson’s ratio, J. Geophys. Res., 70(20), 5249–5257.
  * Walsh (1965b) Walsh, J. B., 1965b.  The effect of cracks on the compressibility of rock, J. Geophys. Res., 70(2), 381–389.
  * Wilson et al. (2003) 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 et al. (2014) 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 (2000) 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 et al. (2007) 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 et al. (1987) 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

Supplementary figures cited in section 4.2.1 : Role of the initial stress
field for plastic deformation.

Figure 13: Figure S1. Plastic deformation for the models with $\Psi=55^{o}$,
S=2, with initial strike-slip stress field (top) and reverse stress field
(bottom). Here C=0. The right panels shows the initial parameters for those
two simulations. Figure 14: Figure S2. Plastic deformation for the models with
$\Psi=15^{o}$, S=1, with initial strike-slip stress field (top) and reverse
stress field (bottom). Here C=0. The right panels shows the initial parameters
for those two simulations. Figure 15: Figure S3. Plastic deformation for the
models with $\Psi=30^{o}$, S=1, with initial strike-slip stress field (top)
and reverse stress field (bottom). Here C=0. The right panels shows the
initial parameters for those two simulations.