# Dynamic evolution of off-fault medium during an earthquake: a micromechanics based model

Marion Y. Thomas, Harsha Bhat

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

## To cite this version:

Marion Y. Thomas, Harsha Bhat. Dynamic evolution of off-fault medium during an earthquake: a micromechanics based model. Geophysical Journal International, 2018, 214 (2), pp.1267-1280. 10.1093/gji/ggy129. hal-03079638

HAL Id: hal-03079638

https://hal.science/hal-03079638

Submitted on 27 Aug 2021

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>
# Dynamic evolution of off-fault medium during an earthquake: a micromechanics based model

Marion Y. Thomas<fnref n="1" /> and Harsha S. Bhat<fnref n="2" />

Received 2018 March 24; in original form 2017 July 29

## SUMMARY

Geophysical observations show a dramatic drop of seismic wave speeds in the shallow off-fault medium following earthquake ruptures. Seismic ruptures generate, or reactivate, damage around faults that alter the constitutive response of the surrounding medium, which in turn modifies the earthquake itself, the seismic radiation and the near-fault ground motion. We present a micromechanics based constitutive model that accounts for dynamic evolution of elastic moduli at high-strain rates. We consider 2-Din-plane models, with a 1-Dright lateral fault featuring slip-weakening friction law. The two scenarios studied here assume uniform initial off-fault damage and an observationally motivated exponential decay of initial damage with fault normal distance. Both scenarios produce dynamic damage that is consistent with geological observations. A small difference in initial damage actively impacts the final damage pattern. The second numerical experiment, in particular, highlights the complex feedback that exists between the evolving medium and the seismic event. We show that there is a unique off-fault damage pattern associated with supershear transition of an earthquake rupture that could be potentially seen as a geological signature of this transition. These scenarios presented here underline the importance of incorporating the complex structure of fault zone systems in dynamic models of earthquakes.

**Key words:** Numerical modelling.

## 1 INTRODUCTION

In the shallow brittle crust, tectonic deformation is believed to occur by seismic or aseismic sliding along active fault planes, controlled by the frictional properties of the fault interface (Scholz 1998). However, faults not only consist of the fine-grained narrow fault core where the extensive shearing is observed, but it is also surrounded by pervasively fractured rocks, within a intricate 3-D geometry (Sibson 1977; Chester et al. 1993; Biegel & Sammis 2004). If fault slip behaviour is intrinsically linked to the properties of the fault interface, the complex structure of fault zone systems impacts the rheological properties of the fault core and the bulk, both of which influence the modes of deformation, and slip, as underlined by recent observations (Andrews 2005; Collettini et al. 2009; Niemeijer et al. 2010; Thomas et al. 2014; Audet & Burgmann 2014). Fault zone structure is, therefore, of key importance to understand the mechanics of faulting.

Micro- and macrostructural field studies, geophysical surveys and laboratory experiments performed on damage zones (e.g. Manighetti et al. 2004; Dor et al. 2006; Savage & Brodsky 2011; Froment et al. 2014; Aben et al. 2017), emphasize these zones as a key component to understand the energy balance of earthquakes

(e.g. Rice 2002; Kanamori 2013). Possible mechanisms responsible for the development of off-fault damage could include the fault geometry, linking of structures, the quasi-static stress field, the process zone associated with fault growth and the coseismic fracture damage (Vermilye & Scholz 1998; Rice et al. 2005; Childs et al. 2009; Faulkner et al. 2011; Vallage et al. 2015). Recent geophysical studies have highlighted the major role played by earthquakes. During seismic faulting, the stress concentration at the tip of the rupture generate, or reactivate, damage (fractures) around faults which can trigger a significant coseismic drops in velocity and modify the constitutive response of the surrounding medium (e.g. Rice et al. 2005). They observe a reduction in elastic stiffness of up to 40%, on spatial scales of hundreds of metres normal to the fault and few kilometres along depth, followed by time-dependent partial recovery over couple of years (Vidale & Li 2003; Hiramatsu et al. 2005; Li et al. 2006; Brenguier et al. 2008; Cochran et al. 2009; Froment et al. 2014). The coseismic change in elastic moduli can influence the rupture dynamics, which has a direct effect on the size of the earthquake, the radiated wavefield and near-fault ground motion (Walsh 1965a,b; Faulkner et al. 2006; Thomas et al. 2017). Moreover, the changes in elastic properties of the off-fault medium controls the amount of elastic strain energy stored during tectonic loading and released during seismic events. The contrast

<footnote n="1">¹Department of Earth Sciences, University of Oxford, South Parks Road, Oxford OX1 3AN, UK. E-mail: mthomas.tectonics@gmail.com</footnote>

<footnote n="2">²Laboratoire de Géologie, Ecole Normale Supérieure, CNRS-UMR8538, PSL Research University, 24 rue Lhomond, Paris, France</footnote>
in elastic moduli, between the fault zones and the host rock, can also induce stress rotation, allowing faults to slip under less optimal far-field stress conditions (Faulkner et al. 2006). Hence, the evolution of damage during earthquake is critical for understanding the nucleation, propagation and arrest of earthquakes.

Several studies have explored the effect of off-fault damage on seismic rupture, using either analytical approaches (Rice et al. 2005; Ngo et al. 2012) or numerical simulations (see following references). The effect of pre-existing damage on the properties of the dynamic rupture (mode, speed and directivity) and final slip has been investigated by prescribing a low-velocity zone around the fault (e.g. Kaneko & Fialko 2011; Cappa et al. 2014; Huang et al. 2014). Another set of models have explored the effect of spontaneous dynamic formation of off-fault damage using a Mohr-Coulomb (e.g. Andrews; 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. 2011; Johri et al. 2014) yield criterion approach, or by modelling off-fault damage as tensile cracks, using a stress-(Yamashita 2000) or fracture-energy-based (Dalguer et al. 2003) criterion. These studies have provided a good insight on the effect a fault zone structure on a dynamic rupture. Nonetheless, they do not include the dynamic changes of elastic properties in the bulk and, therefore, they cannot model the complete feedback between the seismic rupture and off-fault damage. This can be achieved by using an energy-based approach to define the new constitutive law, as developed in the models by Lyakhovsky et al. (1997a), Finzi et al. (2009), Suzuki (2012), Xu et al. (2014) and Lyakhovsky & Ben-Zion (2014). Finally, to reproduce the dynamic evolution of the constitutive response of the bulk, simulations need to duly account for the response of cracks to the applied coseismic loading, using micromechanics-based models. In particular, when focusing on earthquake-related processes, numerical models of off-fault damage need to include a physically motivated crack growth law that takes into account the dependency of loading rate and crack-tip velocities on fracture toughness, as observed in laboratory experiments (Chen et al. 2009; Dai et al. 2010, 2011; Wang et al. 2010, 2011; Zhang & Zhao 2013). The latter constitutes the essential difference between the model used for this study (Bhat et al. 2012; Perol & Bhat 2016; Thomas et al. 2017) and the models aforementioned.

# 2 CONSTITUTIVE MODEL AND NUMERICAL METHOD

## 2.1 General considerations

During an earthquake, ruptures generate/reactivate damage (fractures) around faults whose off-fault fracture density decays with an exponential (Vermilye & Scholz 1998; Wilson et al. 2003; Mitchell & Faulkner 2009; Faulkner et al. 2006, 2011) or power law (Savage & Brodsky 2011) away from the fault. Off-fault damage modifies the microstructure and the constitutive response of the surrounding medium, which in turn, influences the earthquake rupture dynamics, the radiated wavefield and near-fault ground motion. To account for off-fault damage evolution during an earthquake, using a proper micromechanics of fracture, the model must incorporate the following key features. (1) It should allow for accumulation of permanent deformation. (2) It should incorporate a physical criterion for the creation of fresh microcracks and/or the extension of pre-existing ones. (3) It should be consistent with laboratory experiments (e.g. Zhang & Zhao 2013), over a wide range of loading regimes (strain rates ranging from ~$10^{-6}$ s⁻¹ to ~$10^4$ s⁻¹). (4) The presence of

microcracks must trigger a dynamic change in elastic modulus of the bulk with changing stress and strain. The micromechanical formulation used in the model has already been tested and ascertained against laboratory experiments by Bhat et al. (2012), which validates the third point. The following sections provide a brief description of how we incorporate the other features in our constitutive model. Mathematical formulations are provided in the supporting information and an extended, more complete, description of the model can be found in Bhat et al. (2012), Perol & Bhat (2016) and Thomas et al. (2017).

## 2.2 Representation of cracks in the model

In the brittle crust, inelastic deformation occurs by the nucleation, the growth and/or the sliding on pre-existing 'fractures' at different scales. The generic term 'fractures' here encompasses faults and joints but also small-scale features such as microcracks, pores, grain-boundaries, mineral defects, mineral twins, etc. present in all natural rocks. Fractures, under regional tensile loading, can grow when local stresses at the crack tip exceed the rock strength, leading to crack propagation. Locally, an increase in pore pressure can also lead to tensile stresses at the crack tip, and can trigger hydraulic cracking even under compressive loading. Under compressive stress, frictional sliding occurs on pre-existing cracks when the shear stress overcomes the frictional resistance acting on the fracture interface. As the faces slide in opposing direction, they create a tensile wedging force that opens wing cracks at the tips of the shear fracture, which nucleate and grow in $\sigma_1$ direction (most compressive) and open in $\sigma_3$ direction (least compressive).

Following Ashby & Sammis (1990) and Deshpande & Evans (2008), we model the pre-existing flaws present in the medium by penny-shaped cracks of radius $a$ (Fig. 1a). We do not account for nucleation of new cracks. For each numerical simulation, we prescribe a volume density of cracks $N_v$ that remains fixed during the loading (i.e. no nucleation of new cracks), but cracks can grow depending on the local state of stress. For simplicity, only the cracks optimally oriented from a Coulomb friction perspective are taken into account [see Bhat et al. (2011) for a justification]. Therefore, initial flaws are all aligned at the same angle $\Phi = \frac{1}{2} \tan^{-1}(1/f)$ to $\sigma_1$, where $f$ is the coefficient of friction. The density of these initial flaws per unit volume can then be characterized by a scalar $D_0$ defined as:

$$D_0 = \frac{4\pi}{3} N_v (a \cos \Phi)^3 \quad (1)$$

where $a \cos \Phi$ is the projection of the crack radius parallel to the direction of $\sigma_1$.

Then, based on the structural observations described above, we model inelastic deformation during seismic rupture by opening of pre-existing cracks and their propagation (crack growth). We, therefore, define three regimes depending on the overall stress-state: two for compressive loading and one for tensile loading. Hydrofracturing is not yet included in the model. Under compressive loading Regime I, stresses are not high enough to allow sliding along initial flaws. Hence, the solid is assumed to behave like an isotropic linear elastic solid. Regime II is reached when the shear stress $\tau$ overcomes the frictional resistance $f(-\sigma)$ acting on microcracks (normal stress is positive under tensile regime). Inelastic deformation is then accounted for by growing tensile wing cracks at the tip of the penny-shaped cracks (Fig. 1b). If the state of stress turns tensile (Regime III), both penny-shaped cracks and wing cracks can open (Fig. 1b). Hence, in our model, cracks grow parallel to
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Schematics and parameters for simulations of dynamic ruptures in a 2-D in-plane model. We consider a right-lateral planar fault (black line in the middle), embedded in a brittle off-fault medium with a damage-constitutive law. Material properties are defined by the density ($\rho$), the S and P-waves speed ($c_s$ and $c_p$) and the initial damage density ($D_0$). A material contrast across the fault is assumed for $D_0$. Corresponding ‘+’ and ‘-’ directions are defined: ‘+’ highlights the direction of the motion of the more compliant material (Σ1), and ‘-’ of the stiffer material (Σ2). Slip-weakening friction (grey box) acts on a 17$R_0$ or 38$R_0$ km long fault, with $R_0$ corresponding to the process zone size. C and T denote the compressional and the tensional quadrants, respectively. The medium is loaded by uniform background stresses with the maximum compressive stress $\sigma_1$ making an angle of 60° with the fault plane. The thick grey line corresponds to the nucleation prone-patch where the initial shear stress is set up to be just above the fault strength. (b) Mode of inelastic deformation depending on the regime reached during the simulation (see Section2.2).</figcaption></figure>

the $\sigma_1$ axis in the form of wings cracks, each of length $l$ (Fig. 1a). Then the current inelastic state is defined by the scalar $D$ (fraction of volume occupied by microcracks), which describes the current damage state of our solid:

$$D = \frac{4\pi}{3} N_v (a \cos \Phi + l)^3. \qquad (2)$$

Here $D$ approaching 1 corresponds to the coalescence stage that leads to the macroscopic fracture of the solid. Then, to develop the constitutive law, the ‘state’ of the cracks in the medium over time needs to be defined as well. Differentiating eq. (2) with respect to time, leads to the following state evolution law:

$$\frac{dD}{dt} = \left( \frac{3D^{2/3} D_0^{1/3}}{a \cos \Phi} \right) \frac{dl}{dt}, \qquad (3)$$

where $\frac{dl}{dt} \equiv v$ effectively corresponds to the instantaneous wing-crack tip speed.

## 2.3 Physical criterion for dynamic crack growth

At this stage, we now need to define a physical criterion that relates crack growth to the local stress conditions at the tip of the off-fault crack and the ability of the material containing the crack to resist fracture under various loading conditions. In fracture mechanics, this problem is typically solved by requiring that, for a crack to grow, the dynamic microcrack stress intensity factor, $K_I^d$, must overcome the resistance of the material to fracturing, which is given by the dynamic fracture toughness of the material, $K_{IC}$.

Experiments performed on crack initiation, for different rocks types, have highlighted the loading rate dependency of dynamic initiation toughness, $K_{IC}^D$ (Wang et al. 2010, 2011; Zhang & Zhao 2013). Although, little is known about the initiation mechanisms of a crack, based on microstructural observations, it is believed that dislocations, at the defect level, accumulate during loading to subsequently coalesce and form a crack (Sangid 2013). The formation of these dislocation features (or secondary cracks) is controlled by the local stress state, and also how it reaches its critical value, that is, the loading history at the location of the defects (Liu et al. 1998). Once the crack is initiated, the crack growth is controlled by the material properties, that is, the dynamic propagation toughness, $K_{IC}^d$. Experiments on granites, sandstones and marbles have shown it depends on the velocity at which the crack is propagating, that is, $K_{IC}^d(v)$ (Chen et al. 2009; Dai et al. 2010, 2011; Zhou & Aydin 2010; Gao et al. 2015). Both set of experiments (on $K_{IC}^D$ and $K_{IC}^d$) have shown that at high loading rate, when approaching coseismic conditions, it is more difficult to initiate and propagate cracks.

However, most damage models do not properly account for high speed crack growth. They usually use a quasi-static or stress corrosion-like crack growth laws (e.g. Charles Law). In our model, we adopt the crack growth law developed by Bhat et al. (2012) that includes a dynamic criterion for both crack initiation and crack propagation, based on experimental observations. Crack only initiates motion when the stress intensity factor exceeds the dynamic initiation toughness ($K_I^d \ge K_{IC}^D$). Once the crack-tip starts propagating,
the dynamic stress intensity factor, $K_i^d$ must equal the material resistance ahead of the crack-tip for it to grow, that is, $K_i^d = K_{iC}^d(v)$. We, therefore, obtain a nonlinear equation for the crack-tip speed $v$ [eq. (s38) in the supporting information], which in turn is used in eq. (3) to solve for damage evolution. The left-hand side of this equality is derived from the model and depends on the regime, whereas the right-hand side is determined experimentally.

## 2.4 Determining the constitutive stress–strain relationship

To determine the constitutive strain-stress relationship of a damaged solid, we use an energy-based approach, as defined by Rice (1971), Hill & Rice (1973) and Rice (1975). This formalism, thermodynamically argued, determines the inelastic behaviour at macroscopic scale that arises from structural rearrangements at microscale such as the nucleation and growth of microcracks. With this approach, inelastic deformation is treated as a sequence of constrained equilibrium states, defined by the state of stress $\sigma$ and an internal variable, which is here the scalar damage state $D$. Hence, the formalism of equilibrium thermodynamics can be adopted and using the properties of thermodynamic potentials, we create an energetically equivalent solid, where the energy variations takes into account what is happening at the microcrack scale level. Assuming isothermal conditions, the Gibbs free energy density $W^G$ of damaged solid, for a given $\sigma$ and $D$, can then be written as the sum of (1) the free energy $W^{G_e}(\sigma)$ of a solid, without flaws, deforming purely elastically and (2) the free energy $W^{G_i}(\sigma, D)$ corresponding to the contribution of microcracks:

$$W^G(\sigma, D) = W^{G_e}(\sigma) + W^{G_i}(\sigma, D). \quad (4)$$

Then, by fixing $D$ for an imaginary equilibrium state, elastic constitutive laws are applicable and the constitutive relation between these properties and macroscopic stress can be determined in terms of the Gibbs free energy:

$$\epsilon = \frac{\partial W^G(\sigma, D)}{\partial \sigma} = \frac{\partial W^{G_e}(\sigma)}{\partial \sigma} + \frac{\partial W^{G_i}(\sigma, D)}{\partial \sigma} = \epsilon^e + \epsilon^i, \quad (5)$$

where $\epsilon^e$ and $\epsilon^i$ are the associated elastic and inelastic strains, respectively. Hence, with this formalism, local structural rearrangement can be related to corresponding changes in the macroscopic stress and strain states. In practice, internal variables like $D$ will have a time-dependent evolution determined by the local conditions (Section 2.3), but the kinetic aspect of it is taken care of separately (Rice 1971). To compute the additional free energy due to cracks, $W^{G_i}(\sigma, D)$, we need to compute the energy release rate $G$ (crack growth) at the tip of the microcracks, which depends on the stress intensity factor $K_i^d$. As underlined in Section 2.3, $K_i^d$ is stress-dependent and therefore will vary with the regimes. A full description of the constitutive strain-stress relationships is provided by Thomas et al. (2017) and the equations are given in the supporting information.

## 2.5 Numerical method and model setup

We have implemented the constitutive damage model described above in the 2-D spectral element code, SEM2DPACK (Ampuero 2012). A full discussion about the update scheme, the resolution and how we regularize the ‘bimaterial fault’ induced ill-posedness (related to dynamic change of elastic moduli) can be found in Thomas et al. (2017).

To explore the dynamic evolution of off-fault medium during earthquake, we consider a 2-D in-plane model, with a 1-D right lateral fault embedded in a brittle off-fault medium with a damage-constitutive law. To simplify the problem, we assume plane strain conditions. Material properties are defined by the density ($\rho$), the S- and P-waves speed ($c_s$ and $c_p$) and the initial damage density ($D_0$). A material contrast across the fault is assumed for $D_0$. The material on the top part of the fault is always the more compliant for the different studied cases (in grey in Fig. 1a).

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 fault strength $\tau = f(-\sigma^*)$. The parameter $\sigma^*$ corresponds to the modified normal stress as defined by the Prakash & Clifton (1993) law, used to regularize the bimaterial fault ill-posed problem following laboratory experiments (e.g. Wibberley et al. 2008). The friction coefficient $f$, that depends on the cumulated slip $\delta$, drops from a static ($f_s = 0.6$) to a dynamic (= 0.1) value over a characteristic distance $\delta_c$:

$$f = \begin{cases} f_s - (f_s - f_d)\delta/\delta_c & \text{if } \delta \le \delta_c \\ f_d & \text{if } \delta > \delta_c. \end{cases} \quad (6)$$

Slip-weakening friction acts on a $17R_0$ or $28R_0$ long fault, with $R_0$ corresponding to the process zone size. Following Day et al. (2005), the static value of the process zone for a mode II rupture is defined by:

$$R_0 = \frac{9\pi}{32(1 - \nu)} \frac{\delta_c \rho c_s^2}{(f_d - f_s)(-\sigma_{zz})}, \quad (7)$$

where $\nu$ is the Poisson’s ratio, and $\sigma_{zz}$ the stress normal to the fault. Then, to avoid any interference with the propagating dynamic rupture we set the domain to be $5R_0$ perpendicular to the fault, and we applied absorbing boundary conditions on the edge of the computational domain.

For every simulation, we apply uniform background stresses, with the maximum compressive stress $\sigma_1$ making a 60° angle with the fault plane. Initial normal stress ($\sigma_0 = \sigma_{zz}^0$) and shear stress ($\tau_0 = \sigma_{xz}^0$) are uniform along the fault, except for a nucleation-prone patch (thick grey line in Fig. 1) for which $\tau_0$ is set to be slightly above the nominal static strength (~0.03% greater). Reference values for all the different parameters are summarized in Table 1.

## 3 DYNAMIC SIMULATIONS OF OFF-FAULT MEDIUM EVOLUTION

In the following section, we explore the generation of seismically triggered off-fault damage and the associated dynamic evolution of bulk rheology. In particular, we want to investigate the effect of a pre-existing damage zone on the fault rupture propagation and on the evolution of the off-fault medium itself. We start the study with a simple case, a 2-Dright-lateral fault inside an homogeneous medium (Granite), with a small initial damage contrast across the fault. Then, we increase the complexity by introducing an exponential decrease of damage density away from the fault, as described in several field studies (Vermilye & Scholz 1998; Wilson et al. 2003; Mitchell & Faulkner 2009; Faulkner et al. 2006, 2011).

### 3.1 Effect of pre-existing damage on the dynamic evolution of bulk elastic properties

For the first example, we initiate a right-lateral rupture in a damaged medium that presents a contrast in $D_0$ across the fault: $D_{01} = 0.2$ on
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Simulation of a dynamic rupture on a right-lateral fault embedded in a Granite, with a material contrast across the fault in the initial damage density: $D_0 = 0.2$ on the top side of the fault, and $D_0 = 0.1$ on the bottom. (a) Snapshot of the state parameter *D* (density of microcracks in the medium) at *t* = 4.7 s. Damage essentially occurs in the tensional quadrants (*T*) and is more prominent when $D_0$ is higher. Also shown are the '+' and '-' directions (as defined in Fig. 1). (b) Profiles of damage density (red) across the fault with the standard deviation (grey) measured around the location of the profiles (±850 m). Initial value ($D_0$) is display in black. (c) Reduction (in per cent) of *S*-wave and *P*-wave speeds in the medium, in relation to off-fault damage. Formula to compute the new $c_s$ and $c_p$ values, based on the regime, are described in Thomas *et al.* (2017), eqs (37), (38) and (42). For this particular simulation, we record a maximum change of 34.2 per cent for *S*-wave speed and 28.2 per cent for *P*-wave speed.</figcaption></figure>

<table><caption>**Table 1.** Table of parameters for simulations.</caption><thead><tr><td>Parameter</td><td>Symbol</td><td>Value</td></tr></thead><tbody><tr><td>Normal stress on the fault (MPa)</td><td>σ<sub>0</sub></td><td>60.7</td></tr><tr><td>Shear stress on the fault (MPa)</td><td>τ<sub>0</sub></td><td>19.9</td></tr><tr><td>Shear stress on the nucleation patch (MPa)</td><td>τ<sub>n</sub></td><td>36.4</td></tr><tr><td>Static friction coefficient</td><td>f<sub>s</sub></td><td>0.6</td></tr><tr><td>Dynamic friction coefficient</td><td>f<sub>d</sub></td><td>0.1</td></tr><tr><td>Characteristic slip (m)</td><td>δ<sub>c</sub></td><td>1</td></tr><tr><td>Poisson's ratio</td><td>ν</td><td>0.276</td></tr><tr><td>Penny-shaped cracks' radius (m)</td><td>a</td><td>60</td></tr><tr><td>Volume density of cracks (× 10<sup>-7</sup> m<sup>-3</sup>)</td><td>N<sub>v</sub></td><td>1.68 or 3.36</td></tr><tr><td>Density (kg m<sup>-3</sup>)</td><td>ρ</td><td>2.7 × 10<sup>3</sup></td></tr><tr><td>S-wave speed (m s<sup>-1</sup>)</td><td>c<sub>s</sub></td><td>3.12 × 10<sup>3</sup></td></tr><tr><td>P-wave speed (m s<sup>-1</sup>)</td><td>c<sub>p</sub></td><td>5.6 × 10<sup>3</sup></td></tr><tr><td>Process zone size (m)</td><td>R<sub>0</sub></td><td>1054</td></tr></tbody></table>

the top part of the fault (more compliant material), $D_{02} = 0.1$ for the bottom part. The initial flaw size ($a = 60$ m) scales with secondary fractures that usually surround faults that extend over tens of kilometres (e.g. Sowers *et al.* 1994; Vallage *et al.* 2015). This is also a numerical compromise to get off-fault damage for a reasonable computational cost. For an initial microcracks size of 60 m, this leads to a volume density of cracks, $N_v$, of $1.68 × 10^{-7}$ (m$^{-3}$) and $3.36 × 10^{-7}$ (m$^{-3}$), respectively.

Otherwise, the elastic properties of the bulk, $ρ$, $c_s$ and $c_p$, correspond to those of a typical Granite on both sides (Table 1). We apply the slip-weakening friction law over a $17R_0$ long fault.

The final damage density distribution *D* is displayed in Fig. 2(a) for *t* = 4.7 s. The time at which the numerical simulation ends is
chosen to avoid boundary effects. As illustrated by the profiles of damage density (Fig. 2b), a right-lateral seismic rupture along a straight fault triggers a significant rheological contrast across the fault. The response of the damaged elastic solid is different in the compressional and tensional quadrants with more damage in the tensile lobe. We also note that a small difference in the initial damage density significantly impacts the final result. More dynamic damage occurs in the more compliant material ($D_0 = 0.2$). At the maximum, the width of the newly created damage zone reaches $0.9R_0$ and corresponds to the location where the higher slip rate has been recorded (Fig. 4). The extent of the highly damaged zone ($D \ge 0.5$) is about $0.6R_0$ whereas, for the stiffer material, it does not exceed $0.3R_0$. A higher value of $D_0$ implies that the cracks are closer, and therefore interact more with each other as the dynamic rupture propagates along the fault, changing the local state of stress in the medium.

The generation of damage induces a dynamic evolution of the elastic properties which differs from a classic bimaterial rupture. Moreover, the uneven damage density distribution leads to the creation of a spatially heterogeneous medium. Based on eqs (s24), (s25) and (s29) (see supporting information), we record a maximum decrease of 34.2 per cent in S-wave speed and 28.2 per cent in P-wave speed for the top material. In comparison, the changes in S-wave speed is few percent smaller (32.7 per cent) for the stiffer material but not significantly different for the P-wave speed (28.0 per cent). Those results are consistent with geophysical observations of temporal changes in seismic velocity along natural fault following earthquake ruptures. Using seismic and geodetic data from the Calico fault (California), Cochran et al. (2009) have documented a wide damage zone with seismic velocities reduced by 40–50 per cent. Following the 1992 Landers earthquake, Li et al. (1994) have estimated a fault zone width of ~180 m and a strong decrease of shear velocity (~30 per cent). Froment et al. (2014), after the 2008 $M_w$ 6.0 Gofar earthquake (East Pacific Rise), have characterized a damaged fault zone showing 10–20 per cent reduction in P-wave velocity followed by a partial and gradual recovery of the elastic properties within the few months following the earthquakes. Additionally, the data from the seismic survey deployed in Parkfield, 3 months after the 2004 $M_w$ 6.0 main shock, have shown 1.0 to 1.5 per cent decreases in seismic wave velocity within an 200-m-wide zone along the fault strike (Li et al. 2006). However, healing processes took off during the following months, as observed by Brenguier et al. (2008) and Froment et al. (2014). Therefore these values, like Li et al. (2006) suggested, were likely higher right after the earthquake.

## 3.2 Width of the damage zone

The width of the damage zone can provide an observational prediction of the model that could be used to test the model against field data. For the simple case described above, we can see that a small initial difference in damage density actively impacts the final width. On the left tensile quadrant ($D_0 = 0.2$), the newly created damaged zone extends up to $0.9R_0 \approx 950$ m, whereas in the right tensile lobe ($D_0 = 0.1$), it does not exceed $0.5R_0 \approx 530$ m. These values are in the range of field observations. In their study of the fault system in the Shawangunk mountains, Vermilye & Scholz (1998) have recorded a maximum width of 75 m on the extensional part of the fault. Along the Punchbowl fault (San Andreas system), the width of the damage zone on one side was estimated to be $\approx 90$ m (Wilson et al. 2003). In the Atacama fault system, the inferred width of

the damage zone for the tensile quadrant is $\approx 150$ m for the Caleta Coloso fault and $\approx 20$ m for the Blanca fault (Mitchell & Faulkner 2009). Following the 1992 Landers earthquake, Sowers et al. (1994) have mapped mesoscopic fractures that appeared during the seismic event, up to 1.1 km away from the main fault plane. Finally, using optical imagery, Vallage et al. (2015) have inferred a damage zone width of 100–1000 m for the 2013 Balochistan earthquake. Although the prediction of the model is comparable to field data, it is worth underlying few points. First, for some of these studies, the width of the damage zone is likely related to several dynamic events (Vermilye & Scholz 1998; Mitchell & Faulkner 2009; Wilson et al. 2003), whereas the newly observed fractures for the Landers and the Balochistan earthquakes correspond to a single event (Sowers et al. 1994; Vallage et al. 2015). Therefore, the results of the model should be more comparable to the last two cases. Moreover, material contrast (rock type) across natural fault is often observed. Thomas et al. (2017) have shown that seismic events propagating along bimaterial interface generate less damage. For subshear ruptures, the ‘bimaterial effect’ instigates a normal stress change behind the rupture tip, with a tensile stress perturbation in the positive direction and a dynamic compression in the opposite direction (Rice 2002). In turn, the change in normal stress influences the generation of wing cracks, which ultimately reduces the width of the damage zone (by ~20 per cent for this particular case). This should be kept in mind when comparing the observational prediction of the model with field data.

For this model, we can also note the occurrence of ‘branches’, which is a direct consequence of the constitutive law, and does not depend on the resolution as discussed in Thomas et al. (2017). They always appear, even for a better or coarser resolution, and they are related to the initial damage density assigned, as illustrated in the following example (Section 3.6). For this specific case, they make on average a 60° angle with the main fault plane. These type of features have been described in laboratory experiments (Ngo et al. 2012). Nevertheless, caution must prevail, since to accurately capture the localization in numerical simulations, the constitutive laws should have an internal length scale. Because of its complexity, the problem is not currently addressed in the model. Results should, therefore, be taken more qualitatively here, and we are not making any conclusive statement about spacing between branched faults or the width of these localized damage zones.

## 3.3 Crack growth mode

In our model, to generalize the problem for an arbitrary state of stress, as in Drucker & Prager (1952) and Rudnicki & Rice (1975), the normal and shear stresses acting on the crack faces are re-written with their corresponding invariant measures [see Bhat et al. (2012) for a justification]. The normal stress, $\sigma$, is represented by the first invariant of the stress tensor and the shear stress, $\tau$, by the second invariant of the deviatoric stress tensor $S_{ij}$, which gives the distance which the stresses deviate from a state of pure hydrostatic stress:

$$ \sigma = \frac{\sigma_{kk}}{3} \quad \text{and} \quad \tau = \sqrt{\frac{1}{2} S_{ij} S_{ij}} \quad \text{with} \quad S_{ij} = \sigma_{ij} - \sigma \delta_{ij}. \quad (8) $$

Therefore, displaying $\sigma$ and $\tau$ illustrates the stress acting on cracks at one time step (Fig.3). For a right-lateral fault, the rupture tip propagating to the left puts material at the top in tension while the rupture tip on the right induce compression in the medium (T- and C- direction, respectively in Fig.3a). We see the opposite trend
for the bottom material. As a consequence, we observe a consecutive decrease in the first invariant in the upper left and bottom right quadrants. However, the decrease in the first invariant is not sufficient to reach the tensile **Regime III** ($\sigma - \sigma_0 / |\sigma_0| > 1$), and damage is not occurring through crack opening (Fig. 3d). We observe a decrease in the second invariant essentially in the compressional quadrant (Fig. 3b). However, it is actually the combined effect of an increase in $\sigma$ (negative in compression) and a reduction in $\tau$ that triggers damage. **Regime II** is reached when the shear stress $\tau$ overcomes the frictional resistance $f(-\sigma)$ acting on microcracks (Fig. 3d). The inelastic deformation is then accounted for by

growing tensile wing cracks at the tip of the penny-shaped cracks (cf. Fig. 1b).

## 3.4 Influence of dynamically evolving bulk on seismic rupture

Fig. 4 displays the cumulative slip, slip rate and normalized stress on the fault, with a time increment of 0.7 s, for the model presented above. We compare this simulated case (coloured lines) with a right-lateral rupture occurring in a pure elastic medium presenting the same properties apart from the initial damage (thin black curves).

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Normalized 1st (a) and 2nd (b) stress invariants, $\sigma$ and $\tau$, at $t = 4.1$ s, for a dynamic rupture on a right-lateral fault embedded in a Granite, with a material contrast across the fault for $D_0$ (cf. Fig. 2). Although we observe a decrease in the first invariant, it does not reach the tensile **Regime III** ($\sigma - \sigma_0 / |\sigma_0| > 1$). (c) gives the regime under which damage was generated. **Regime II** is reached when the shear stress $\tau$ overcomes the frictional resistance $f(-\sigma)$ acting on microcracks. Therefore, under these stress field conditions, damage occurs by growing tensile wing cracks. Slip rate on the fault (black/white curves) is superimposed on the snapshots.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Simulation of a dynamic rupture on a right-lateral fault embedded in a Granite, with a material contrast across the fault for $D_0$ (cf. Fig. 2). Cumulative slip (a), normalized normal stress change (b) and slip rate (c) on the fault are displayed with a time increment of 0.7 s. Coloured curves correspond to the dynamic simulation with the damage evolution law, dotted grey curves represent a simulation with the same parametrization but for a pure elastic medium.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Logarithm of the absolute value of fault-normal particle velocities $v_n$, at $t = 4.1$ s, for a dynamic rupture on a right-lateral fault, with (a) off-fault damage evolution and (b) for a pure elastic medium. Both simulations share the same parametrization, but for the initial damage contrast across the fault for the model (a), as described in Fig. 2. (c) Displays the differences between the two models. Plotting $log(|v_n|)$ emphasized the high frequency content related to the propagation of seismic waves in the damaged off-fault medium. Corresponding fault slip rate (black/red curves) is superimposed on the snapshots.</figcaption></figure>

Looking at the cumulative slip (Fig. 4a), little difference is noted between the pure elastic mode and the simulation with dynamic damage. In both cases, rupture propagates as a crack, and only a small decrease in cumulative slip is observed at the rupture front for each time steps. We can also observe that the cumulative slip in the negative direction is slightly smaller than in the positive direction, where less damage is recorded. Fig. 4(c) displays the normal stress on the fault for different time steps. We observe positive and negative undulations behind the rupture front, which indicate tensile and compressive changes relative to the background compressive stress state, respectively. This is because damage generation locally changes the state of stress and therefore has an impact on both shear and normal traction acting on the fault.

Evolution of slip rate along the fault plane strongly differs from a classic elastic model (Fig. 4b). Although the rupture is bilateral in both the cases, we observe the development of slip rate oscillations for the simulation with off-fault damage. We also note that these oscillations increase in amplitude as the width of the damage zone becomes larger. As shown by the spacing between symbols, which correspond to the value at each node, oscillations in slip rate are well resolved numerically. They are likely related to the dynamic evolution of elastic properties behind the rupture front that creates a low velocity zone (LVZ). Indeed, in their study of the 1992 Landers earthquake, Li et al. (1994) have recognized the presence of trapped waves, in relation to the presence of a ~180 m wide damage fault zone, where a strong decrease of shear velocity (~30 per cent) was observed. The material contrast (LVZ) can produce internal wave reflections which, in turn, can influence rupture propagation. However, with regard to the complex pattern of the LVZ created in

our model, it is difficult to decipher between the different parameters that can influence the intricate feedbacks we observe (e.g. velocity contrast, width and relative distance between branches, etc). Finally, it is important to note that we observe almost no modulations of the rupture front compared to the simulation with a pure elastic medium. The likely explanation for this simple damage case is that the dynamic rupture, which propagates at subshear velocity on average (~2.7 km s⁻¹), interacts with an homogeneous material. However, the radiated waves may interact with the LVZ behind the rupture front and can further interfere with the rupture front itself in some cases.

## 3.5 Impact of damage on particle velocity field

Figs 5(a) and (b) display the logarithm of the absolute value of fault normal velocity $v_n$, at $t = 4.1$ s, for both the pure elastic case and the simulation with off-fault damage generation. Fig. 5(c) displays the differences between the two models. Plotting $log(|v_n|)$ emphasizes the high frequency content related to the propagation of seismic waves in the damaged off-fault medium in comparison to the elastic case. In particular, we observe that the high frequency signal essentially arises behind the rupture front when the slip rate starts to decrease. However, if we look at the difference between the two models (Fig.5c), we can also note undulations in the normal velocity field ahead of the rupture, which are likely related to the damage generation triggered by the seismic waves. Indeed, a comparison with Fig. 3(c) emphasizes the parabolic features we can observe at both edges of the simulated box (around −8$R_0$ and 8$R_0$).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Simulation of a dynamic rupture on a right-lateral fault embedded in a Granite, with an exponential decrease in from $D_0=0.5$ to 0.1 on the top side of the fault, and $D_0 = 0.1$ on the bottom. (a) Snapshot of the state parameter *D* at *t* = 6.5 s. (b) Corresponding profiles across the fault. *D* is displayed in red, the standard deviation (±1.8 km) in grey and the initial value ($D_0$) in black. Damage essentially occurs in the tensional quadrants, but also in the compressional quadrant on the top side of the fault (where $D_0$ is higher), essentially because of seismic waves propagating in the medium (see Fig. 7). (c) and (d) gives the cumulative slip and slip rate on the fault are with a time increment of 0.7s. Coloured curves correspond to the dynamic simulation with the damage evolution law. The dotted grey lines represent a simulation with the same parametrization but for a pure elastic medium.</figcaption></figure>

The complexity in fault slip rate, related to the off-fault damage generation, are likely responsible for the high frequency content of the fault-normal particle velocities. This is combined with the effect of seismic waves propagating in a off-fault damage medium, where elastic properties have been modified dynamically. These observations are consistent with laboratory experiments observations, where (Passelegue et al. 2016) have related the high frequency radiation observed during laboratory earthquakes with the amount of damage produced. Similar features are also observed on near-fault strong motion records of real earthquakes (Housner 1947; Wald & Heaton 1994; Semmane et al. 2005; Dunham et al. 2011). However, many factors can contribute to the high-frequency content: fault-roughness and/or slip heterogeneities, which also result in ground acceleration spectra that are flat at high frequency (Dunham et al. 2011). An upcoming paper by the authors will explore the combined effects of off-fault damage and fault roughness on the radiated ground motion.

## 3.6 Complex damage zone: effect of an exponential decrease of $D_0$ on seismic rupture under subshear and supershear regime

### 3.6.1 Geometry of the newly created damage zone

For the second numerical experiment, we increase the complexity of our model by introducing an exponential decrease of damage density away from the fault, as described in several field studies (Vermilye & Scholz 1998; Wilson et al. 2003; Mitchell & Faulkner 2009; Faulkner et al. 2006, 2011). On the top side of the fault, the

initial damage density varies from $D_0=0.5$ to $D_0=0.1$ over a distance equivalent to the process zone $R_0 \sim 1$ km. The values chosen for $D_0$ do not precisely reflect the measured values in the field as more observations are needed. The goal of this relatively simple scenario is to evaluate qualitatively the influence of a pre-existing damage zone, accumulated over several seismic events, on the following earthquake rupture and on the dynamic evolution of the bulk. However, the width of the initial damage zone and the initial crack size (60 m) are constrained by the field observations following the 1992 Landers and 2013 Balochistan earthquakes (Sowers et al. 1994; Vallage et al. 2015). They both described mesoscopic fractures (10 to ~150 m long) up to 1 km away from the main fault plane. For comparison, we assigned a uniform value of initial damage density ($D_0=0.1$) for the bottom part of the model. Compared to the previous example, we also extend the length of the fault to $38R_0$.

The damage density distribution *D* at the end of the numerical simulation (*t* = 6.8 s) is displayed in Figs 6(a) and (b). The pattern of the initial damage density highly influences the development of the damage zone. The dynamic rupture generates branches in the tensional quadrant of the stiffer material whereas the softer material displays an increase in *D* in both, the positive and negative direction. Nevertheless, the rupture traveling on the tensional side activates and/or interacts more with the off-fault damage, creating a homogeneous pattern, with branches only appearing at the edges of the highly damaged zone. This is likely related to the higher initial value of damage density in the softer material ($D_0=0.5$ near the fault), as observed for the previous case (Section 3.1). Moreover, the profiles of damage density (Fig. 6b) show that the width of the zone where $D=1$ seems to reach a steady-state between $-7R_0$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Temporal evolution of a dynamic rupture occurring on a right-lateral fault embedded in a Granite, with an exponential decrease in initial damage density from $D_0 = 0.5$ to $D_0 = 0.1$ over $R_0$ on the top side of the fault, and $D_0 = 0.1$ on the bottom (see Fig. 6 for the final stage). (a) Regime under which damage is generated at several time steps. (b) Corresponding snapshots of the norm of the particles velocity ($\sqrt{v_n^2 + v_p^2}$). Fault slip rate (red curves) is superimposed for every time steps. The transition to the supershear regime [$S = 1.19 < 1.77$] does not occur at the same time on both sides (at $t \sim 4.8$ s on the negative side, at $t \sim 5.4$ s on the positive side). It has an effect on the shape of the damage zone (see also Fig. 6a). We observe the development of a Mach cone on the negative side at $t \sim 6.5$ s. However, due to the complexity of the dynamically evolving bulk, it occurs slightly behind the rupture front, and does not appear yet on the positive side.</figcaption></figure>

and $-11.5R_0$. On the contrary, the width of the LVZ constantly increases in the stiffer material as the crack-like rupture grows in subshear regime (Fig. 6d). Then, because of the parameters assigned for the slip-weakening friction law, the rupture goes from subshear to supershear (Fig. 6d). The parameter $S = (\tau_p - \tau_0)/(\tau_0 - \tau_r)$, defined by Andrews (1976), where $\tau_p = -f_s\sigma_n$ and $\tau_r = -f_d\sigma_n$ are, respectively, the peak and residual strength, gives the threshold at which the seismic rupture can become supershear ($S< 1.77$). In our simulation, $S = 1.19$ and the transition to the supershear regime highly impacts the damage generation on both sides.

Fig. 7 displays the temporal evolution of damage and particle velocity in the medium. As the rupture propagates below the Rayleigh wave speed, damage is generated behind the rupture front in the tensional quadrants, with more damage where $D_0$ is higher. Tensile wing cracks are also opening ahead of the rupture front, where the S-wave field concentrates (Fig. 7a), although it leads to a much smaller D than behind the rupture front (Fig. 6a). The transition to the supershear regime, which occurs at different time on both sides, has a direct impact on the shape of the damage zone. On the negative side, a new pulse is generated ahead of the rupture at $x \approx 10.5R_0$, and at $x \approx 12.5R_0$ on the positive side. In both the cases, it

coincides with a decrease in the width of the damage zone (Fig. 6a). This is likely related to the decrease in the stress intensity factor as the rupture speed increases. The stress intensity factor of the main rupture in mode II, $K_{II}^d$, for a dynamically growing crack ($v > 0$), is expressed as the product of a crack-tip velocity dependent function times the value of the static stress intensity factor for a crack of a given length $l$, $K_{II}$:

$$K_{II}^d(v) \approx \frac{1 - v/c_R}{\sqrt{1 - v/c_p}} K_{II} \quad \text{with} \quad K_{II} = \varphi(\sigma, \tau)\sqrt{\pi L}, \quad (9)$$

where $c_R$ corresponds to the Rayleigh wave speed. The function depends on the given applied loading ($\sigma, \tau$). As a consequence, when approaching the Rayleigh wave speed, $K_{II}^d$ becomes small (even if the crack length has increased quite a bit), resulting in modest off-fault stress perturbation and hence little change in damage. The location where the width of the LVZ is the smallest in the softer material ($x = 13R_0$) closely corresponds to the location where two distinct pulses are being generated (Figs 6d and 7a), when the rupture becomes supershear. As shown by the snapshot at $t = 5.6$ s, when the first pulse is not strong enough to change optimally the state of stress, damage only occurs behind the front of the pulse propagating below the Rayleigh wave speed. As in the previous
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Logarithmic damage density distribution as a function of perpendicular distance from the fault. (a) and (b) display the upper left and bottom right quadrant for the simulation of a dynamic rupture with a material contrast across the fault in $D_0$: $D_0 = 0.2$ on the top side of the fault, and $D_0 = 0.1$ on the bottom. See Fig. 2(a) for the location of the profiles. (c) and (d) display the upper left and upper right quadrant for the simulation of a dynamic rupture with an exponential decrease in $D_0$ from 0.5 to 0.1 over $R_0$ on the top side of the fault, and $D_0 = 0.1$ on the bottom. See Fig. 6(a) for the location of the profiles. The grey boxes give the initial damage density distribution, $D_0$, at the beginning of the simulations.</figcaption></figure>

time steps, we can also observe damage fields on both sides related to the concentration of the *S*-waves. Moreover, on the negative side, wing cracks are opening ahead of the newly created pulse. It is more diffuse on the positive side because the transition to the supershear regime has not yet happened and the highest value of $D_0$ is in the compressive quadrant.

At $t = 6.5$ s, both rupture fronts have developed a second pulse ahead. However, looking at the particle velocity (Fig. 7b) we only observe the development of a Mach cone on the negative side. This is related to the complexity of the dynamically evolving bulk. With respect to the local changes in elastic moduli, only the left front went supershear before the simulation ended. Concurrently, damage is generated at three locations on the negative side: at the rupture front propagating at the *S*-wave speed, and behind the first and second pulses in the tensional quadrant. This corresponds to an increase in the width of the LVZ (Fig. 6a). On the positive side, damage occurs behind the rupture front propagating below the Rayleigh wave speed and because of seismic waves propagating in the softer medium.

### 3.6.2 Impact of the complex damage zone on the seismic rupture

Figs 6(c) and (d) display, respectively, the cumulative slip and slip rate on the fault, with a time increment of 0.7 s. We compare the model (coloured lines) with a right-lateral rupture occurring in a pure elastic medium (dotted grey lines). In terms of cumulative slip, the fault has slipped slightly less than for the elastic case when more damage has accumulated. However, the differences are minimal. On the other hand, evolution of slip rate along the fault plane strongly differs from a classic elastic model and, more importantly, complex $D_0$ distribution leads to a bilateral but asymmetric rupture (Fig. 6d). Compared to the elastic case, damage evolution in the bulk counterbalances the systematic increase of slip rate related to the slip weakening law. Then the rupture propagating on the positive side, because of the smaller $D_0$ in the tensional quadrant, activates/interacts less with the cracks in the off-fault medium. As a consequence, in spite of the observed slip rate oscillations related to the development of an LVZ, we note little modulation of the rupture front, compared to the elastic case. On the contrary, on the negative side, because of the complex interaction between the evolving off-fault medium and the dynamic event, the rupture front propagates
at a slower rate than the elastic simulation until it reaches the super-
shear regime. This transition also happens earlier for the simulation
with off-fault damage because of the reduction in $c_s$ speed in the
medium.

4 DISCUSSION AND CONCLUSION

In this paper, we have numerically investigated the role of seismic
events on the dynamic evolution of the off-fault medium and un-
derlined the intricate feedback that affects concomitantly the earth-
quake rupture processes. These features cannot be reproduced with
models allowing for elastoplastic deformation. Moreover, the con-
stitutive law developed in this model, allows to properly account for
the microphysics of damage evolution related to earthquake rupture,
by including the dependency of fracture toughness loading rate and
crack-tip velocities (Chen et al. 2009; Dai et al. 2010, 2011; Wang
et al. 2010, 2011; Zhang & Zhao 2013). This constitutes the essen-
tial difference with the pre-existing models that seek to reproduce
the dynamic evolution of the bulk (e.g. Lyakhovsky et al. 1997b;
Xu et al. 2014).

The two scenarios presented here underline the importance of incorporating the complex structure of fault zone systems in dynamic models of earthquakes. A small difference in initial damage actively impacts the final damage pattern. The second numerical experiment, in particular, highlights the complex feedback that exists between the evolving medium and the seismic event. An interesting test would be to use optical image correlation for past events, such as the 1999 Izimit or the 2002 Denali earthquakes, to see if we can observe a signature of the transition to the supershear regime in the off-fault damage generation. Applying the method used for the 2013 M<sub>w</sub> 7.7 Balochistan earthquake (Vallage et al. 2015), one can use the fault parallel and fault normal displacements to look at the along strike variations of the damage zone width.

As a comparison with geological survey of damage fault zones,
Fig. 8 displays the logarithmic damage density distribution as a
function of fault normal distance, before and after a dynamic event.
For the first example, we initiate a right-lateral rupture in a simple
medium that presents a contrast in $D_0$ across the fault: $D_{01} = 0.2$
on the top part of the fault, $D_{02} = 0.1$ for the bottom part. After just
one dynamic event, we can see that to the first order, when $D_{01} =
0.2$ (Fig. 8a) the generated damage distribution is compatible with
the field observation of an exponential decrease of crack density
away from the fault (Vermilye & Scholz 1998; Wilson et al. 2003;
Mitchell & Faulkner 2009; Faulkner et al. 2006, 2011). For a smaller
value of $D_0$ (Fig. 8b), the distribution is closer to a power law
decrease, as observed by Savage & Brodsky (2011). A higher value
of $D_0$ could represent a fault zone that has undergone more seismic
events.

For the second numerical experiment, we increase the complexity
of our model by introducing an pre-existing exponential decrease of
damage density away from the fault. The first diagram (Fig.8c) sam-
ples the fault zone in the tensile quadrant. Although we started with
an exponential decrease in $D_0$, some of the final profiles deviate
from the initial distribution. However, our model does not take into
account other mechanisms that can kick off during the interseismic
period and heal, at least partly, the newly created damage. Indeed,
geophysical observations suggest that the damage effect is transient,
with gradual (sometimes incomplete) recovery of the elastic proper-
ties (e.g. Brenguier et al. 2008; Froment et al. 2014). This evolution
is likely related to healing processes that affect microcracks, frac-
tures and faults through precipitation of soluble materials or clay

mineralization (Mitchell & Faulkner 2008). Therefore, one should
expect the cracks to heal faster if they are more connected (higher
value of *D*). In the compressional quadrant, however (Fig.8c), even
with an increase in damage density, the generated damage distri-
bution still displays an exponential decrease of crack density away
from the fault. Although, it is worth noticing that for the profiles
that record the highest increase in *D*, we tend to lose the ‘stair-case’
distribution of the initially damage state.

Finally, to be consistent with the geophysical observations that
suggest a gradual recovery of the elastic properties (e.g. Froment
et al. 2014), one should develop models that take into account the
complex competition between the intensity of the coseismic rupture,
the efficiency of healing processes and the recurrence time between
earthquakes. A significant step further, very challenging however,
would then be to develop a constitutive law that accounts for the
evolution of elastic properties inside the fault zone due to both the
dynamic damage and ‘healing/sealing’ between dynamic events.

ACKNOWLEDGEMENTS

This study was supported by the Agence National de la Recherche
(ANR) GeoSMEC contract ANR-12-BS06-0016. The paper was
completed at Oxford with support for MT from the Natural Envi-
ronment Research Council large grant, Looking Inside the Conti-
nents from Space (NE/K011006/1). We are grateful to Eric Dun-
ham and an anonymous reviewer for their comments that improved
our manuscript. We also thank J.-P. Ampuero for providing the
SEM2DPACK code (available at https://sourceforge.net/projects/s
em2d/) and for his useful comments on how to add the new consti-
tutive model to the SEM code.

REFERENCES

Kanamori, H. & Rivera, L., 2013. *Earthquakes: Radiated Energy and the Physics of Faulting*. 170, Abercrombie, R., McGarr, A., Di, Toro G. & Kanamori, H., 3–13, *Energy Partitioning During an Earthquake*.

Aben, F.M., Doan, M.-L., Gratier, J.-P. & Renard, F., 2017. Coseismic damage generation and pulverization in fault zones, in *Fault Zone Dynamic Processes*, Thomas, M.Y., Mitchell, T.M. & Bhat, H.S., pp. 47–80, eds John Wiley & Sons, Inc.

Ampuero, J.P., 2012. *Sem2dpack, a spectral element software for 2d seismic wave propagation and earthquake source dynamics, v2.3.8 [data set], zenodo*, http://doi.org/10.5281/zenodo.230363, Tech. rep., Caltech.

Andrews, D.J., 1976. Rupture velocity of plane strain shear cracks, *J. geophys. Res.*, **81**(32), 5679–5687.

Andrews, D.J., 2005. Rupture dynamics with energy loss outside the slip zone, *J. geophys. Res.*, **110**(B1), B01307, doi.org/10.1029/2004JB003191.

Ashby, M.F. & Sammis, C.G., 1990. The damage mechanics of brittle solids in compression, *Pure appl. Geophys.*, **133**(3), 489–521.

Audet, P. & Burgmann, R., 2014. Possible control of subduction zone slow-earthquake periodicity by silica enrichment, *Nature*, **510**(7505), 389–392.

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

Bhat, H.S., Rosakis, A.J. & Sammis, C.G., 2012. A micromechanics based constitutive model for brittle failure at high strain rates, *J. Appl. Mech. Trans. ASME*, **79**(3), 031016, doi:10.1115/1.4005897.

Bhat, H.S., Sammis, C.G. & Rosakis, A.J., 2011. The micromechanics of westerley granite at large compressive loads, *Pure appl. Geophys.*, **168**(12), 2181–2198.

Biegel, R.L. & Sammis, C.G., 2004. Relating fault mechanics to fault zone structure, in *Advances in Geophysics*, Vol. 47, pp. 65–111, Elsevier.
Brenguier, F., Campillo, M., Hadzioannou, 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.

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

Chen, R., Xia, K., Dai, F., Lu, F. & Luo, S., 2009. Determination of dynamic fracture parameters using a semi-circular bend technique in split Hopkinson pressure bar testing, *Eng. Fract. Mech.*, **76**(9), 1268–1276.

Chester, F.M., Evans, J.P. & Biegel, R.L., 1993. Internal structure and weakening mechanisms of the San Andreas fault, *J. geophys. Res.*, **98**(B1), 771–786.

Childs, C., Manzocchi, T., Walsh, J.J., Bonson, C.G., Nicol, A. & Schöpfer, M.P., 2009. A geometric model of fault zone and fault rock thickness variations, *J. Struct. Geol.*, **31**(2), 117–127.

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., Viti, C., Smith, S.A.F. & Holdsworth, R.E., 2009. Development of interconnected talc networks and weakening of continental low-angle normal faults, *Geology*, **37**(6), 567–570.

Dai, F., Chen, R., Iqbal, M. & Xia, K., 2010. Dynamic cracked Chevron notched Brazilian disc method for measuring rock fracture parameters, *Int. J. Rock Mech. Min. Sci.*, **47**(4), 606–613.

Dai, F., Xia, K., Zheng, H. & Wang, Y., 2011. Determination of dynamic rock mode-i fracture parameters using cracked Chevron notched semi-circular bend specimen, *Eng. Fract. Mech.*, **78**(15), 2633–2644.

Dalguer, L.A., Irikura, K. & Riera, J.D., 2003. Simulation of tensile crack generation by three-dimensional dynamic shear rupture propagation during an earthquake, *J. geophys. Res.*, **108**(B3), doi.org/10.1029/2001JB001738.

Day, S.M., Dalguer, L.A., Lapusta, N. & Liu, Y., 2005. Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture, *J. geophys. Res.*, **110**(B12), B12307, doi:10.1029/2005JB003813.

Deshpande, V.S. & Evans, A.G., 2008. Inelastic deformation and energy dissipation in ceramics: a mechanism-based constitutive model, *J. Mech. Phys. Solids*, **56**(10), 3077–3100.

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 appl. Geophys.*, **163**(2), 301–349.

Drucker, D.C. & Prager, W., 1952. Soil mechanics and plastic analysis or limit design, *Q. Appl. Math.*, **10**(2), 157–165.

Dunham, E.M., Belanger, D., Cong, L. & Kozdon, J.E., 2011. Earthquake ruptures with strongly rate-weakening friction and off-fault plasticity, part 2: Nonplanar faults, *Bull. seism. Soc. Am.*, **101**(5), 2308–2322.

Faulkner, D.R., Mitchell, T.M., Healy, D. & Heap, M.J., 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), doi.org/10.1029/2010JB007788.

Finzi, Y., Hearn, E.H., Ben-Zion, Y. & Lyakhovsky, V., 2009. Structural properties and deformation patterns of evolving strike-slip faults: numerical simulations incorporating damage rheology, *Pure appl. Geophys.*, **166**(10), 1537–1573.

Froment, B., McGuire, J.J., van der Hilst, R.D., Gouedard, P., Roland, E.C., Zhang, H. & Collins, J.A., 2014. Imaging along-strike variations in mechanical properties of the gofar transform fault, east Pacific rise, *J. geophys. Res.*, **119**(9), 7175–7194.

Gabriel, A.A., Ampuero, J.P., Dalguer, L.A. & Mai, P.M., 2013. Source properties of dynamic rupture pulses with off-fault plasticity, *J. geophys. Res.*, **118**(8), 4117–4126.

Gao, G., Yao, W., Xia, K. & Li, Z., 2015. Investigation of the rate dependence of fracture propagation in rocks using digital image correlation (DIC) method, *Eng. Fract. Mech.*, **138**(Complete), 146–155.

Hill, R. & Rice, J.R., 1973. Elastic potentials and structure of inelastic constitutive laws, *SIAM J. Appl. Math.*, **25**(3), 448–461.

Hiramatsu, Y., Honma, H., Saiga, A., Furumoto, M. & Ooida, T., 2005. Seismological evidence on characteristic time of crack healing in the shallow crust, *Geophys. Res. Lett.*, **32**(9), doi.org/10.1029/2005GL022657.

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, *Geophys. Res. Lett.*, **37**, L02306, doi.org/10.1029/2009GL041888.

Housner, G.W., 1947. Characteristics of strong-motion earthquakes, *Bull. seism. Soc. Am.*, **37**(1), 19–31.

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

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, *J. geophys. Res.*, **119**(2), 1251–1272.

Kaneko, Y. & Fialko, Y., 2011. Shallow slip deficit due to large strike-slip earthquakes in dynamic rupture simulations with elasto-plastic off-fault response, *Geophys. J. Int.*, **186**(3), 1389–1403.

Liu, C.S., Liu, S.Y., Lallemand, S.E., Lundberg, N. & Reed, D.L., 1998. Digital elevation model offshore Taiwan and its tectonic implications, *Terre Atmos. Ocean. Sci.*, **9**(4), 705–738.

Li, Y.-G., Aki, K., Adams, D., Hasemi, A. & Lee, W.H.K., 1994. Seismic guided waves trapped in the fault zone of the landers, California, earthquake of 1992, *J. geophys. Res.*, **99**(B6), 11705–11722.

Li, Y.-G., Chen, P., Cochran, E.S., Vidale, J.E. & Burdette, T., 2006. Seismic evidence for rock damage and healing on the San Andreas fault associated with the 2004 m 6.0 Parkfield earthquake, *Bull. seism. Soc. Am.*, **96**(4B), S349–S363.

Lyakhovsky, V. & Ben-Zion, Y., 2014. A continuum damage-breakage faulting model and solid-granular transitions, *Pure appl. Geophys.*, **171**(11), 3099–3123.

Lyakhovsky, V., Ben-Zion, Y. & Agnon, A., 1997a. Distributed damage, faulting, and friction, *J. geophys. Res.*, **102**(B12), 27635–27649.

Lyakhovsky, V., Reches, Z., Weinberger, R. & Scott, T.E., 1997b. Non-linear elastic behaviour of damaged rocks, *Geophys. J. Int.*, **130**(1), 157–166.

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

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

Mitchell, T.M. & Faulkner, D.R., 2008. Experimental measurements of permeability evolution during triaxial compression of initially intact crystalline rocks and implications for fluid flow in fault zones, *J. geophys. Res.*, **113**(B11), B11412.

Mitchell, T.M. & Faulkner, D.R., 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, *J. Struct. Geol.*, **31**(8), 802–816.

Ngo, D., Huang, Y., Rosakis, A., Griffith, W.A. & Pollard, D., 2012. Off-fault tensile cracks: a link between geological fault observations, lab experiments, and dynamic rupture models, *J. geophys. Res.*, **117**(B1), doi.org/10.1029/2011JB008577.

Niemeijer, A., Marone, C. & Elsworth, D., 2010. Fabric induced weakness of tectonic faults, *Geophys. Res. Lett.*, **37**, L03304, doi.org/10.1029/2009GL041689.
Palmer, A.C. & Rice, J.R., 1973. The growth of slip surfaces in the pro-
gressive failure of over-consolidated clay, *Proc. R. Soc. A*, **332**(1591),
527–548.

Passelègue, F.X., Spagnuolo, E., Violay, M., Nielsen, S., Di Toro, G. & Schubnel, A. 2016. Frictional evolution, acoustic emissions activity, and off-fault damage in simulated faults sheared at seismic slip rates, *Journal of Geophysical Research-solid Earth*, **121**(10), 7490–7513.

Perol, T. & Bhat, H.S., 2016. Micromechanics-based permeability evolution in brittle materials at high strain rates, *Pure appl. Geophys.*, **173**(8), 2857–2868.

Prakash, V. & Clifton, R., 1993. Time resolved dynamic friction measure-
ments in pressure shear, in *Experimental Techniques in the Dynamics of
Deformable Solids*, Vol. **165**, pp. 33–48.

Rice, J.R., 1971. Inelastic constitutive relations for solids - an internal-
variable theory and its application to metal plasticity, *J. Mech. Phys.
Solids*, **19**(6), 433–455.

Rice, J.R., 1975. Continuum mechanics and thermodynamics of plasticity in relation to microscale deformation mechanisms, in *Constitutive Equations in Plasticity*, chap.2, ed. Argon, A., MIT Press.

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, *Bull. seism. Soc. Am.*, **95**(1), 109–134.

Rudnicki, J.W. & Rice, J.R., 1975. Conditions for localization of deformation in pressure-sensitive dilatant materials, *J Mech Phys Solids*, **23**(6), 371–394.

Sangid, M.D., 2013. The physics of fatigue crack initiation, *Int. J. Fatigue*, **57**, 58–72.

Savage, H.M. & Brodsky, E.E., 2011, Collateral damage: evolution with displacement of fracture distribution and secondary fault strands in fault damage zones, *J. geophys. Res.*, **116**, doi:10.1029/2010JB007665.

Scholz, C.H., 1998. Earthquakes and friction laws., *Nature*, **391**, 37–42.

Semmane, F., Cotton, F. & Campillo, M., 2005. The 2000 tot-
tori earthquake: A shallow earthquake with no surface rupture
and slip properties controlled by depth, *J. geophys. Res.*, **110**(B3),
doi.org/10.1029/2004JB003194.

Sibson, R.H., 1977. Fault rocks and fault mechanisms, *J. Geol. Soc. Lond.*, **133**(3), 191–213.

Sowers, J.M., Unruh, J.R., Lettis, W.R. & Rubin, T.D., 1994. Relationship of the kickapoo fault to the Johnson valley and Homestead valley faults, San Bernardino county, California, *Bull. seism. Soc. Am.*, **84**(3), 528–536.

Suzuki, T., 2012, Understanding of dynamic earthquake slip behavior using damage as a tensor variable: Microcrack distribution, orientation, and mode and secondary faulting, *J. geophys. Res.*, **117**(B5), doi.org/10.1029/2011JB008908.

Templeton, E.L. & Rice, J.R., 2008. Off-fault plasticity and earthquake rupture dynamics: 1. dry materials or neglect of fluid pressure changes, *J. geophys. Res.*, **113**(B9), doi.org/10.1029/2007JB005529.

Thomas, M.Y., Avouac, J.-P., Gratier, J.-P. & Lee, J.-C., 2014. Lithological control on the deformation mechanism and the mode of fault slip on the longitudinal valley fault, Taiwan, *Tectonophysics*, **632**, 48–63.

Thomas, M.Y., Bhat, H.S. & Klinger, Y., 2017. Effect of brittle off-fault damage on earthquake rupture dynamics, in *Fault Zone Dynamic Processes: Evolution of Fault Properties During Seismic Rupture*, Vol. **227**, pp. 255–280, eds Thomas, M.Y., Mitchell, T.M. & Bhat, H.S., John Wiley & Sons, Inc.

Vallage, A., Klinger, Y., Grandin, R., Bhat, H.S. & Pierrot-Deseilligny, M., 2015. Inelastic surface deformation during the 2013 mw 7.7 Balochistan, Pakistan, earthquake, *Geology*, **43**(12), 1079–1082.

Vermilye, J.M. & Scholz, C.H., 1998. The process zone: A microstructural view of fault growth, *J. geophys. Res.*, **103**(B6), 12223–12237.

Vidale, J.E. & Li, Y.G., 2003. Damage to the shallow landers fault from the nearby Hector mine earthquake, *Nature*, **421**(6922), 524–526.

Wald, D.J. & Heaton, T.H., 1994. Spatial and temporal distribution of slip for the 1992 landers, California, earthquake, *Bull. seism. Soc. Am.*, **84**(3), 668–691.

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.

Wang, L., Hainzl, S., Sinan, Özeren M. & Ben-Zion, Y., 2010. Postseismic deformation induced by brittle rock damage of aftershocks, *J. geophys. Res.*, **115**(B10), B10422, doi:10.1029/2010JB007532.

Wang, Q., Feng, F., Ni, M. & Gou, X., 2011. Measurement of mode i and mode ii rock dynamic fracture toughness with cracked straight through flattened Brazilian disc impacted by split Hopkinson pressure bar, *Eng. Fract. Mech.*, **78**(12), 2455–2469.

Wibberley, C.A., Yielding, G. & Di Toro, G., 2008. Recent advances in the understanding of fault zone internal structure: A review, *Geol. Soc. Lond. Spec. Publ.*, **299**(1), 5–33.

Wilson, J.E., Chester, J.S. & Chester, F.M., 2003. Microfracture analysis of fault growth and wear processes, Punchbowl fault, San Andreas system, California, *J. Struct. Geol.*, **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 appl. Geophys.*, **172**(5), 1243–1267.

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

Zhang, Q. & Zhao, J., 2013. Effect of loading rate on fracture toughness and failure micromechanisms in marble, *Eng. Fract. Mech.*, **102**, 288–309.

Zhou, X. & Aydin, A., 2010. Mechanics of pressure solution seam growth and evolution, *J. geophys. Res.*, **115**(B12), B12207, doi.org/10.1029/2010JB007614.

SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Table S1.** Parameters of the damage constitutive model.

**Figure S1.** Synthetic seismograms of fault-parallel (a) and fault-normal (b) velocities with the corresponding Fourier amplitude spectra (FAS). Black curves correspond to the dynamic simulation with homogeneous elastic properties but different initial damage (case 1, Fig. 2). Coloured curves correspond to a simulation with the same parametrization but within a pure elastic medium. ‘x’ and ‘y’ coordinates give the location of the synthetic seismometers. Modified from Thomas *et al.* (2017).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for this paper.