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

S U M M A R Y Geophysical observations show a dramatic drop of seismic wave speeds in the shallow offfault 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-D in-plane models, with a 1-D right 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.


I N T RO D U C T I O N
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 2006). 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 per cent, 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 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 & Harris 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. 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 . 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. 2010Dai et al. , 2011Wang et al. 2010Wang et al. , 2011Zhang & 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.

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. 2006Faulkner et al. , 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 −1 to ∼10 4 s −1 ). (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).

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 σ 1 direction (most compressive) and open in σ 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 = 1 2 tan −1 (1/ f ) to σ 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: where acos is the projection of the crack radius parallel to the direction of σ 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 τ overcomes the frictional resistance f(−σ ) 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 the σ 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: 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: where dl/dt ≡ v, effectively corresponds to the instantaneous wingcrack tip speed.

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 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 d I , must overcome the resistance of the material to fracturing, which is given by the dynamic fracture toughness of the material, K D I C . Experiments performed on crack initiation, for different rocks types, have highlighted the loading rate dependency of dynamic initiation toughness, K D I C (Wang et al. 2010Zhang & 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 d I C . Experiments on granites, sandstones and marbles have shown it depends on the velocity at which the crack is propagating, that is, K d I C (v) (Chen et al. 2009;Dai et al. 2010Dai et al. , 2011Zhou & Aydin 2010;Gao et al. 2015). Both set of experiments (on K D I C and K d I C ) 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 d I ≥ K D I C ). Once the crack-tip starts propagating, the dynamic stress intensity factor, K d I must equal the material resistance ahead of the crack-tip for it to grow, that is, K d I = K d I C (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.

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 σ 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 σ and D, can then be written as the sum of (1) the free energy W Ge (σ ) of a solid, without flaws, deforming purely elastically and (2) the free energy W G i (σ , D) corresponding to the contribution of microcracks: 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: where e and 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 (σ , 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 d I . As underlined in Section 2.3, K d I 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.

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 damageconstitutive law. To simplify the problem, we assume plane strain conditions. Material properties are defined by the density (ρ), the Sand 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 (Fig. 1a).
Rupture propagation along the fault plane is governed by a slipweakening friction law (e.g. Palmer . Slip occurs when the on-fault shear stress reaches the fault strength τ = f(−σ *). The parameter σ * 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 δ, drops from a static (f s = 0.6) to a dynamic (f d = 0.1) value over a characteristic distance δ c : 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: where ν is the Poisson's ratio, and σ 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 σ 1 making a 60 • angle with the fault plane. Initial normal stress (σ 0 = σ 0 zz ) and shear stress (τ 0 = σ 0 xz ) are uniform along the fault, except for a nucleationprone patch (thick grey line in Fig. 1) for which τ 0 is set to be slightly above the nominal static strength (∼0.03 per cent greater). Reference values for all the different parameters are summarized in Table 1.

DY N A M I C S I M U L AT I O N S O F O F F -FAU LT M E D I U M E V O L U T I O N
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-D right-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. 2006Faulkner et al. , 2011.

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 0 1 = 0.1 on the top part of the fault (more compliant material), D 0 2 = 0.2 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 ≥ 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 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 6.0 main shock, have shown 1.0 to 1.5 per cent decreases in seismic wave velocity within an 200-mwide 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.

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 950 m, whereas in the right tensile lobe (D 0 = 0.1), it does not exceed 0.5R 0 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 90 m (Wilson et al. 2003). In the Atacama fault system, the inferred width of the damage zone for the tensile quadrant is 150 m for the Caleta Coloso fault and 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 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.

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, σ , is represented by the first invariant of the stress tensor and the shear stress, τ , 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: Therefore, displaying σ and τ 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  Fig. 2). Although we observe a decrease in the first invariant, it does not reach the tensile Regime III ((σ − σ 0 )/ σ 0 > 1). (c) gives the regime under which damage was generated. Regime II is reached when the shear stress τ overcomes the frictional resistance f(−σ ) 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. 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 (σ − σ 0 )/ σ 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 σ (negative in compression) and a reduction in τ that triggers damage. Regime II is reached when the shear stress τ overcomes the frictional resistance f(−σ ) 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). 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 rightlateral rupture occurring in a pure elastic medium presenting the same properties apart from the initial damage (thin black curves).

Influence of dynamically evolving bulk on seismic rupture
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 −1 ), 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.

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 −8R 0 and 8R 0 ). 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 ?) 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.

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. 2006Faulkner et al. , 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 ∼ 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 an 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 and −11.5R 0 . On the contrary, the width of the LVZ constantly  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. 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 = (τ p − τ 0 )/(τ 0 − τ r ), defined by Andrews (1976), where τ p = −f s σ n and τ r = −f d σ 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 ≈ 10.5R 0 , and at x ≈ 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 d I I , 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 : where c R corresponds to the Rayleigh wave speed. The function depends on the given applied loading (σ , τ ). As a consequence, when approaching the Rayleigh wave speed, K d I I 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 snaphsot 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 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.

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 offfault medium and the dynamic event, the rupture front propagates at a slower rate than the elastic simulation until it reaches the supershear regime. This transition also happens earlier for the simulation with off-fault damage because of the reduction in c s speed in the medium.

D I S C U S S I O N A N D C O N C L U S I O N
In this paper, we have numerically investigated the role of seismic events on the dynamic evolution of the off-fault medium and underlined the intricate feedback that affects concomitantly the earthquake rupture processes. These features cannot be reproduced with models allowing for elastoplastic deformation. Moreover, the constitutive 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. 2010Dai et al. , 2011Wang et al. 2010Wang et al. , 2011Zhang & Zhao 2013). This constitutes the essential 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 w 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 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 = ∼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. 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 0 1 = 0.1 on the top part of the fault, D 0 2 = 0.2 for the bottom part. After just one dynamic event, we can see that to the first order, when D 0 = 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. 2006Faulkner et al. , 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) samples 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 properties (e.g. Brenguier et al. 2008;Froment et al. 2014). This evolution is likely related to healing processes that affect microcracks, fractures 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 distribution 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.

A C K N O W L E D G E M E N T S
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 Environment Research Council large grant, Looking Inside the Continents from Space (NE/K011006/1). We are grateful to Eric Dunham 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 constitutive model to the SEM code.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online. Table S1. Parameters of the damage constitutive model. Figure S1. Synthetic seismograms of fault-parallel (a) and faultnormal (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.