Off-fault Plasticity and Earthquake Rupture Dynamics: 1. Dry Materials or Neglect of Fluid Pressure Changes

[ 1 ] We analyze inelastic off-fault response during earthquakes. Spontaneous crack-like rupture, with slip weakening, is modeled in 2-D plane strain using an explicit dynamic finite element procedure. A Mohr-Coulomb type elastic-plastic description describes the material bordering the fault. We identify the factors which control the extent and distribution of off-fault plasticity during dynamic rupture. Those include the angle with the fault of the maximum compressive prestress, the seismic S ratio, and the closeness of the initial stress state to Mohr-Coulomb failure. Plastic response can significantly alter the rupture propagation velocity, delaying or even preventing a transition to supershear rupture in some cases. Plastic straining also alters the residual stress field left near the fault. In part 1, we consider ‘‘dry’’ materials bordering the fault, or at least neglect pore pressure changes within them. Part 2 addresses the effects of fluid saturation, showing that analysis procedures of this part can describe undrained fluid-saturated response. Elastic-plastic laws of the type used are prone to shear localization, resulting in an inherent grid dependence in some numerical solutions. Nevertheless, we show that in the problems addressed, the overall sizes of plastic regions and the dynamics of rupture propagation seem little different from what are obtained when we increase the assumed plastic hardening modulus or dilatancy parameter above the theoretical threshold for localization, obtaining a locally smooth numerical solution at the grid scale. Evidence for scaling of some localization features with a real (nongrid) length scale in the model is also presented.


Observations
[2] Geologic observations of the structure surrounding mature, large-displacement faults show that the typical fault zone consists of an inner core composed of finely granulated material bordered by a gouge layer that grades into fractured-damaged wall rock [Chester and Logan, 1986;Chester et al., 1993;Caine et al., 1996;Biegel and Sammis, 2004;Chester et al., 2004]. In many faults this typical structure, as shown in Figure 1, is not symmetric about the fault [Chester et al., 2004], as might especially be the case when the principal fault surface separates different host rocks [Wibberley and Shimamoto, 2003;Dor et al., 2006]. Dor et al. [2006] and Shi and Ben-Zion [2006] suggest that asymmetries in damage are associated with variations in lithology of the wall rock, in that a contrast in lithology across a fault may cause dynamic rupture events to have a preferred direction.

Previous Modeling
[3] The elastically predicted stresses surrounding a rapidly propagating rupture can increase significantly over the surrounding prestress, particularly for shear rupture approaching its ''limiting speed'' (the Rayleigh wave speed for mode II rupture and the shear wave speed for mode III). Poliakov et al. [2002] and Rice et al. [2005] examined elastically predicted stress fields near propagating ruptures to determine where they violated the Mohr-Coulomb failure criterion and thereby roughly estimated patterns of off-fault damage activation. They found that the zones of potential damage activation depend strongly on the initial prestress orientation, rupture propagation direction, and propagation speed. Rice et al. [2005] also considered pore fluid effects in the form of undrained poroelastic material response and use of Terzaghi effective stress in the failure criterion. Figure 2 shows the 2-D geometry of the shear ruptures considered. Those studies showed that Y, the angle of most compressive initial principal stress to the fault, controls the location with respect to the fault of potential inelastic deformation, with steep Y favoring inelasticity on the extensional side and shallow Y favoring inelasticity on the compressional side. Also, pore fluid effects enhanced failure on the compressional side while reducing it on the extensional side of the fault. The extent of those potential zones of failure increase as rupture velocity, v r , increases and becomes comparable to the size of the low-speed, low-stress drop, slip-weakening zone length, R 0 as v r approaches c R , the Rayleigh wave speed. For the set of seven earthquakes studied by Heaton [1990], Rice et al. [2005] estimate average values of R 0 of 20 to 40 m at midseismogenic depth, with a range of 1 to 80 m.
[4] The high elastically predicted stresses near the propagating rupture were thus found inconsistent with the assumption that off-fault response is elastic for a brittle material with, presumably, little or no cohesion, that is typical of the damaged regions bordering fault zones. Yamashita [2000], Dalguer et al. [2003], Ando and Yamashita [2007], Andrews [2005], Ben-Zion and Shi [2005], and Duan [2008] have investigated dynamic rupture propagation with spontaneous generation or activation of damage in the off-fault material. Yamashita [2000], using the finite difference method, and Dalguer et al. [2003], using the discrete element method, modeled the generation of off-fault damage as the formation of tensile cracks. They considered stress states with Y = 45°a nd found that zones of tensile cracks occur along the extensional side of the fault and increase in extent with increasing rupture propagation distance along the main fault. The boundary integral method was used by Ando and Yamashita [2007] to model the spontaneous dynamic formation of mesoscopic-scale fault structures during earthquake rupture. Using the finite difference method and incorporating a nondilatant Mohr-Coulomb type yield criterion, Andrews [2005] confirmed that for Y = 45°, plastic flow in predamaged (i.e., cohesionless) material accumulates along the extensional side of the fault, and the fault normal extent of plastic strain is proportional to the propagation distance, at least for the distances that it was feasible to study.

Objectives of the Present Work
[5] We investigate factors controlling the theoretical extent, magnitude, and distribution of off-fault inelastic deformation, and the resulting effects of that deformation on the dynamics of shear rupture, by extending the elastic analyses of Poliakov et al. [2002] and Rice et al. [2005] to full elastic-plastic analysis of crack-like rupture, like in the work by Andrews [2005]. We extend the work of Andrews [2005] by including effects of the prestress direction, in the form of Y, and allowing also for the possibility of strain hardening and plastic dilatancy. We also identify spontaneous strain localization in off-fault plastic regions which has gone unnoticed in some previous work. In this first paper, part 1, we consider the case of ''dry'' materials bordering the fault, or at least we neglect any mechanical effects of dynamic pore pressure alterations. In part 2, Viesca et al. [2008] discuss the effects of fluid saturation on the phenomena addressed here. They show that analyses of rupture propagation with undrained elastic-plastic deformation off the fault can largely be described by the same elastic-plastic formulation that we use here, but with altered values assigned for the elastic bulk modulus and plastic constitutive parameters.
[6] Propagating shear rupture is analyzed in our work using the dynamic finite element method in the form of ABAQUS/Explicit [ABAQUS, Inc., 2005]. We first address what level of cohesion in the off-fault yield condition is needed to suppress plastic deformation during rupture. We then investigate, for material with zero cohesion, the factors which control the location, extent, and magnitude of offfault inelastic deformation. Those factors are the seismic S ratio, which measures proximity of the initial shear stress along the fault to peak frictional strength there, the angle Y of the most compressive prestress, and the closeness CF of the initial stress state to Mohr-Coulomb off-fault failure. Next we discuss how the creation of plastic strain during dynamic rupture affects the evolution of rupture velocity, the residual level of fault-parallel compressive stress, and the amount of slip on the fault. We also analyze the localizations of plastic strain that sometimes occur, confirming that the parameter range of plastic properties allowing their occurrence is consistent with the existing theoretical background, and we suggest how such features can be controlled in future analyses. Although the strain localizations observed have an inherent grid size dependence, we show that their  Model geometry, sense of slip, and initial stresses (which result in an angle Y of the most compressive stress to the fault) for numerical studies of shear rupture. Four quadrants, two compressional (minuses) and two extensional (pluses), are defined by the first motions experienced in the material surrounding the fault. The terms compressional and extensional also refer to the signs of the fault-parallel strain xx along the fault walls near the rupture front.
presence does not significantly alter the overall size of the plastic region, or dynamics of rupture propagation, from cases for which the choice of modestly different values of plastic hardening and dilatancy precludes localized plastic deformation.

Constitutive Response of Fault and Adjoining Material
[7] We aim to characterize the inelastic deformation occurring in ''dry'' material surrounding a shear crack propagating along a fault in a 2-D configuration, as in Figure 2. The initial stress state is given by with normal stresses positive in tension. The fault plane is parallel to, and the slip and propagation directions perpendicular to, the intermediate principal stress direction of the tectonic prestress, s ij 0 . If pore fluid is present, we regard these as initial effective stresses. In part 1, the material is assumed to, effectively, have no pore fluids, but Viesca et al. [2008] consider the effects of pore fluid pressure changes in part 2, by transforming all present constitutive parameters into new parameters by rules based on the work by Rudnicki [1984aRudnicki [ , 1984bRudnicki [ , 2000].

Slip-Weakening Friction
[8] The shear strength of a fault can be approximately described using a slip-weakening Coulomb friction law, where the strength is proportional to the normal stress on the fault, s n , and weakens with increasing slip, Du, on the fault. Proposed by Ida [1972] and Palmer and Rice [1973], slip-weakening friction is widely used as a failure criterion to describe the earthquake rupture process, although it does not include possibly severe rate weakening at seismic slip rates or fully describe other dynamic weakening processes considered in current research on fault zone physics. In the simplest case of linear degradation of the friction coefficient with slip Du, the shear strength t of the fault is where t p and t r are the peak and residual shear strength, respectively, given as t p = f s (Às n ) and t r = f d (Às n ) with static and dynamic coefficients of friction f s and f d , respectively. D c is the critical slip-weakening distance over which the degradation of strength occurs. The length of the static slip-weakening zone at incipient failure can be estimated, using the expression of Palmer and Rice [1973], as when the Poisson ratio is n = 0.25; G is the elastic shear modulus and G is the fracture energy (G = (t p À t r ) D c /2 for linear slip weakening).
[9] Equation (3) is valid in the limit of large S (i.e., s xy 0 only slightly greater than t r ), where S is defined as Andrews [1976] observed in numerical simulations that when the seismic S ratio is sufficiently low, S < 1.77, the rupture propagation speed can transition from the sub-Rayleigh to the supershear regime by formation of a daughter crack ahead of the main rupture. Rice [1980] showed that in a large S range, the ratio of the actual slipweakening zone size R to R o is a universal function of rupture velocity which decreases monotonically to zero at the Rayleigh wave speed. Analogous results for general S are given by Rice et al. [2005] and in the supershear range of rupture propagation by Bhat et al. [2007a].
[10] Linear slip weakening is adopted here as a simple description for the fault. However, the appropriateness of a linear form has been questioned in recent observational studies [Abercrombie and Rice, 2005], and lab-and theorybased studies of friction at seismic slip rates emphasize an inherent rate dependence that is not represented in the slipweakening concept [Rice, 2006]. Abercrombie and Rice [2005] confirmed the result of Ohnaka [1996] that linear slip weakening can match observationally inferred earthquake source parameters over a wide range of magnitudes only when D c depends on the final slip along the fault, a problematic result because it implies that the early weakening behavior of the fault depends on the final slip that it will sustain.

Elastic-Plastic Off-Fault Material Description
[11] Materials such as rocks and soils exhibit pressuredependent yielding [Brace et al., 1966;Mogi, 1972Mogi, , 1974Jaeger and Cook, 1976;Hirth and Tullis, 1992] in which the onset of plastic deformation depends on the mean normal stress at temperatures too low and/or timescales too short for creep. Inelastic deformation in brittle rocks under compressive stress occurs primarily as frictional sliding on fissure surfaces and microcracking. We follow Rudnicki and Rice [1975] in formulating a pressure-dependent elastic-plastic constitutive relation to describe the general features observed in the deformation of granular materials and brittle rocks. The constitutive description we use can be understood in terms of the behavior of a material subjected to a hydrostatic compressive stress, Às (with s positive in tension), and a shear stress, t. The relationships between stress and strain increments for a material with linear elastic isotropy can be written as during elastic-plastic response, where G, K and h are the shear, bulk, and plastic hardening moduli, respectively. A plastic strain increment is d p g = dg À dt/G. The dilatancy factor, b, is the ratio of an increment in volumetric plastic strain to an increment in shear plastic strain: d p = bd p g.
The plastic hardening modulus is defined such that hd p g = dt + mds when d p g > 0 (d p g < 0 is disallowed). During plastic response, the ratio of (dt + mds) to h is always positive; for hardening both (dt + mds) and h are positive, for softening both are negative, and for ideally plastic response both vanish such that their limit exists as a positive quantity corresponding to d p g.
[12] Rudnicki and Rice [1975] generalize the material description for arbitrary stress states by taking s = s kk /3 and t = t ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1=2Þs ij s ij p , the second invariant of the deviatoric stress s ij = s ij À d ij s kk /3. The yield condition proposed by Drucker and Prager [1952], a pressuredependent modification of the Huber-von Mises yield criterion, is a simple choice for describing granular or cracked materials. The material description above reduces to it when m is taken to be independent of stress level s, so that the yield surface defining the shear strength as a function of mean normal stress is a straight line with slope m. The Drucker-Prager yield criterion [Drucker and Prager, 1952;Lubliner, 1990, chapter 3.3

.3] is given by
Here b is the cohesion, or shear strength at zero mean normal stress. A more detailed description of the elasticplastic formulation, including the treatment of plastic dilational and shear straining, can be found in Appendix A.
[13] In plane strain, the Drucker-Prager (DP) yield criterion approximates the Mohr-Coulomb (MC) criterion ( where t and s n are the shear and normal traction on any plane, tan f is the internal friction coefficient and c is the cohesion. The DP and MC criteria coincide exactly for 2-D stress states like those in equation (1) when the out of plane principal stress is given by s zz = (s xx + s yy ) /2 = s kk /3 (i.e., when s zz = 0). For those stress states, the cohesion and friction coefficients are related by b = c cos f and m = sin f, which we take to define b and m here. Both the MC and DP yield criteria are idealizations of the behavior of brittle materials subjected to compressive stresses. According to available multiaxial stress experiments, neither is fully suitable for describing the inelastic behavior of rocks [Davis and Selvadurai, 2002;Colmenares and Zoback, 2002]. We use the DP description, which is widely used and well tested for geomechanics applications, for this study since it has the implementation advantage over the MC description that the yield surface is smooth. The DP model is an available material description in the explicit finite element package that we use, although through the user subroutine VUMAT, descriptions for other elastic-plastic models, including MC, can be introduced.

Numerical Model
[14] We consider a 2-D mode II configuration in which the rupture propagation direction and slip direction coincide and are in the x direction, with the fault along the x axis ( Figure 2). A 2-D rectangular mesh is used, composed of four-noded reduced integration plane strain elements (type CPE4R in ABAQUS) with linear relations between strains and displacement gradients. The reduced integration procedure uses a one-point integration, based on the standard uniform strain formulation [Flanagan and Belytschko, 1981], to form the element stiffness matrix. This formulation, in which the local stress state s ij assigned to a given element is based on the strain history at a single (centroidal) integration point, is adopted because, in addition to providing a significant computational advantage over a full integration scheme, for elastic-plastic materials, elements that are more precisely integrated tend to propagate artificial constraints [Nagtegaal et al., 1974]. Hourglass controls based on the artificial stiffness method derived by Flanagan and Belytschko [1981] are used within ABAQUS/Explicit to prevent the reduced integration procedure from leading to zero-energy modes, or ''hourglassing,'' which, if allowed, would dominate the solution. A set of overlapping nodes In 2-D plane strain, with s zz = 1 2 (s xx + s yy ), the criteria are equivalent, and the friction angles and cohesions can be related by the simple equations shown above. The closeness of the stress state to MC or DP failure is defined as CF, and along stress paths that maintain s zz = s kk /3, which holds initially, CF is the same for MC and DP.
(split nodes) is present along the fault. Along the predefined fault, we use a split node contact procedure (ABAQUS user subroutine VFRIC) to prescribe the weakened shear strength resulting from slip during rupture along the fault. During slip, tangential forces are applied at each node along the fault, consistent with the shear strength of the fault (equation (2)). Details of the implementation are provided in Appendix B.
[15] Initial stresses in the material surrounding the fault are prescribed to be consistent with an assumed angle Y (Figure 2), and the out of plane principal stress is s zz 0 = s kk 0 /3, in which case the MC and DP failure criteria coincide (along paths with s zz = s kk /3). In the analyses and presentation of results that follow, all stresses have a nondimensionalization given by the initial fault-normal stress, s yy 0 , and lengths have a nondimensionalization R 0 (equation (3)). The element width Dx is chosen so that the static slip-weakening zone size is well resolved, with Dx = R 0 /20 (i.e., typically 1 to 2 m in real physical units), and sometimes much smaller.
[16] Along a portion of the fault of length L 0 c , we impose an initial shear stress distribution consistent with a slipweakening shear crack in energetic equilibrium like in the work by Kame et al. [2003] and Bhat et al. [2004] to nucleate a dynamic shear rupture. The length of the nucleation zone, L 0 c , is chosen somewhat greater than the critical nucleation length at instability, calculated by Palmer and Rice [1973] and given here for n = 0.25 as for the large S limit coinciding with singular elastic crack mechanics with small-scale yielding. This initial alteration in shear stress initiates dynamic rupture at both ends of the nucleation zone at the start of the simulation to produce a bilateral right-lateral shear rupture along the predefined fault.
[17] Absorbing elements (called ''infinite elements'' in ABAQUS) surround the entire mesh to minimize reflections from the boundaries. These elements introduce normal and shear tractions on the boundary of the finite element mesh that are proportional to the normal and shear components of velocity at the boundary, with damping constants chosen as the wave impedance factors to minimize reflections of dilational and shear wave energy. Forces are applied between the boundary of the plane strain elements and the infinite elements consistent with the prescribed initial stress state.
[18] We focus most of our investigation on nondilatant (b = 0) plastic response with no hardening (h = 0) of the yield surface. We do, however, investigate briefly the effects of nonzero b and h, in part to resolve issues of strain localization [Rudnicki and Rice, 1975] which occur with the h = 0 and b = 0 model (and would occur with modestly

Results and Discussion
[19] We address the following questions: Where does offfault inelastic deformation develop and how does that deformation affect the dynamics of earthquake rupture? We investigate the role of the initial stress state, in the form of Y, the seismic S ratio, and the closeness CF of the initial stress state to the elastic-plastic yield criterion, in determining the location and extent of off-fault plastic strain during a dynamic shear rupture. Y together with S and the static and dynamic friction coefficients, f s and f d , along the fault completely determine the 2-D initial stress state used (equation (1)), normalized by the initial fault-normal stress, (Às yy 0 ), as follows: The parameter, CF, is Figure 6. The equivalent plastic shear strain field during dynamic shear rupture is plotted for Y ranging from 10°to 45°. Nucleation lengths, L 0 c , are identical for each case, but the times at which the plots are made vary. The fault friction coefficients are f s = 0.5 and f d = 0.1, and the off-fault elastic-plastic material properties are tan f = 0.6, c = 0, b = 0, and h = 0. The grid size is Dx = R 0 /20. The location of equivalent plastic shear strain during rupture propagation depends strongly on Y, with inelastic deformation occurring primarily on the compressional side for Y 20°, and exclusively on the extensional side for Y ! 45°. Note that for fixed S, f s , and f d , varying Y leads to a variation in CF.
is the radius of the Mohr circle for stress, s m = (s xx + s yy )/2 is the mean stress, and tan f and c are the internal friction coefficient and cohesion, respectively, as defined for the MC yield criterion (along stress paths that maintain s zz = s kk /3, which holds initially, the value of CF is the same for DP). Figure 3 shows a schematic of the initial stress state and its distance from the yield surface along a line of constant mean normal stress (Às kk /3) for the MC and DP yield criteria. If CF > 1 the initial state of stress violates the yield criterion. When the fault friction coefficients and off-fault MC material properties are specified, CF depends only on S and Y, and the three cannot be varied independently.
[20] In the results that follow, we first describe how S, Y, and CF control the plastic strain field that develops during dynamic shear rupture. Next we discuss how the resulting stress field, rupture velocity, and slip accumulation are different from what would be predicted during a shear rupture with off-fault elastic response. In sections 4.1 -4.5, the resulting plastic strain during dynamic rupture are presented as contours of ''equivalent'' plastic strain, g eq p :  (see equation (A5)). Figure 4 shows the evolution of shear stress and equivalent plastic strain during a crack-like bilateral dynamic shear rupture. Plastic strain develops behind the rupture, with the active zone of plastic strain, d p g/dt > 0, located immediately behind the rupture tip.

Amount of Cohesion Necessary to Prevent Plastic Strain
[21] We first investigate how much cohesion is necessary to prevent plastic deformation altogether in the material surrounding the fault. Using the finite element method, we calculate the dynamic change in stress, s ij À s ij 0 , surrounding a mode II shear rupture propagating in the sub-Rayleigh regime as it accelerates toward the Rayleigh speed at a time when the rupture velocity, v r , reaches 0.86c s . Because the stress concentration ahead of a mode II crack depends on the rupture velocity, the calculated level of cohesion required is sufficient to prevent plastic yielding only for sub-Rayleigh ruptures with v r 0.86c s . No fixed finite cohesion, c, can prevent off-fault yielding if v r is allowed to get arbitrarily close to the Rayleigh speed [Rice, 1980], in the sense that finite slip then accumulates over a vanishingly small slipweakening zone R. We use that change in stress to calculate the total stress field surrounding a dynamic shear rupture, for Y ranging from 5°to 85°and S ranging from 0.5 to 2.5, for rupture along a fault with static and dynamic coefficients of friction of f s = 0.6 and f d = 0.1. We determine the amount of cohesion, c, required to prevent violation of the MC yield criterion For each Y À S combination, the minimum c necessary to satisfy the inequality in equation (13) is calculated for tan f = 0.6. The amount of cohesion, c/(Às yy 0 f s ), needed to prevent plastic deformation has a strong Y dependence ( Figure 5) and is lowest for Y between 10°and 25°. Because the level of cohesion required depends on the fault-normal stress, at greater depths in the Earth more cohesion will be necessary to prevent off-fault plasticity since normal stress increases with depth. If the cohesion of the material surrounding the fault is independent of depth over some interval, this means, paradoxically, that plastic deformation in the material surrounding a fault will occur more readily as depth increases in that interval. Nevertheless, if yield occurs at all depths, the plastically activated zone may be larger at shallower depths if D c has weak or no depth dependence . In sections 4.2 -4.5, we assume that the material surrounding the fault is highly cracked or granulated from the fault evolution process and from stressing in previous earthquakes so that it has effectively lost all cohesion. That neglects cohesion from mineralization processes during the interseismic period.

Location of Plastic Strain With Respect to the Fault
[22] The first motions experienced in the material surrounding a fault due to P wave arrivals define four quadrants, two compressional and two extensional, as shown in Figure 2. We find that the angle of most compressive stress determines the locations, with respect to the fault, where plastic strain will occur, as anticipated from Poliakov et al. [2002], Kame et al. [2003], and Rice et al. [2005]. Plastic deformation occurs primarily on the compressional side of the fault for shallow angles of most compressive stress and on the extensional side for higher angles of most compressive stress. Figure 6 demonstrates the effect of Y, for Y ranging from 10°to 45°with the seismic S ratio fixed at S = 1.0 and off-fault elastic-plastic material properties, tan f = 0.6 and c = 0, on the location of plastic strain with respect to the fault with contours of ''equivalent'' plastic shear strain. As Y increases from 10°to 20°, the plastic deformation shifts from occurring exclusively on the compressional side of the fault to occurring with equal extent on the compressional and extensional sides. As Y increases further, the location of plastic strain shifts to the extensional side of the fault, with all plastic deformation occurring on the extensional side for Y ! 45°. During rupture, the zone of active deformation occurs near the rupture tip. In particular, for the cases shown in Figures 7a and 7b, the instantaneous dg p /dt is nonzero along the sharp line between the zero and nonzero g eq p ahead of the rupture tip and the sharp line between the intermediate g eq p ( 0.9t p /(2G)) and high g eq p (!1.7t p /(2G)) immediately behind the rupture tip.

Extent of Plastic Zone
[23] The seismic S ratio and the closeness CF of the offfault material to failure control the extent of the zone of plastic deformation, for a given rupture propagation distance, and the magnitude of equivalent plastic strain within that zone. The extent of the zone of plastic deformation increases with increasing CF and decreasing S. Equation (11), which expresses CF in terms of the initial stress state, can be rewritten for c = 0 as S, Y, and CF cannot be chosen independently of one another. CF increases with decreasing S (i.e., with increased closeness to failure on the fault itself) for a given Y. Figure 8 shows the relationship between CF, S, and Y for f s = 0.6 and f d = 0.1 on the fault, and off-fault material with tan f = 0.6 (the fault and surrounding material have the same peak frictional strength), which are consistent with values    [Brudy et al., 1997;Zoback and Harjes, 1997].
[24] The extent of plastic deformation increases with decreasing seismic S ratio, as shown in Figures 7 and 9 for Y = 14°(with f s = 0.6, f d = 0.1, and tan f = 0.6) and Y = 56°(with f s = 0.45, f d = 0.045, and tan f = 0.6), respectively. As S decreases, the extent of the plastic zone increases due, in part, to the material surrounding the fault being closer to failure initially, with CF increasing from 0.79 to 0.93 for Y = 14°and from 0.47 to 0.65 for Y = 56°, as S decreases from 1.5 to 0.75.
[25] Strong evidence of strain localization is observed for Y = 56°in Figures 9a and 9b. These features remain present with increasing grid refinement and their spacing, at least at the grid resolution shown, is set by the grid size. Rudnicki and Rice [1975] found that materials with pressure-dependent yielding can have localization of deformation (which corresponds dynamically to planar deformation waves with zero propagation speed) even during quasi-static deformation with (sufficiently small) positive hardening. We will discuss the localization features we observe subsequently.

Effects on Rupture Dynamics
[26] Simulations of crack-like dynamic shear rupture in a homogenous elastically deforming material show that, far behind the leading edge of a propagating rupture, the change in fault-parallel stress s xx is a small fraction of the stress drop, with s xx becoming more compressive on the compressive side of the fault, and more tensile on the extensional side of the fault. In our analyses of dynamic rupture in material with elastic-plastic response, we instead observe a significant residual s xx within the zone of plastic deformation left behind the rupture front. Figure 10 shows s xx becomes less compressive in a zone of plastic deformation on the compressional side of the fault while it becomes more compressive in a zone of plastic deformation on the extensional side, irrespective of the angle of most compressive stress. The differences between the residual s xx on the compressional and extensional sides of the fault, if they could be measured, could be used as an indicator of rupture directivity in the most recent event since fault-parallel stresses in the zone of inelastic deformation will be less compressive on the compressional side of the fault and more compressive on the extensional side. This is at least so if s xx is equal on both sides before rupture; that need not be the case. A complete study of this point, not provided here, will require attention to the effects of multiple prior ruptures, including cases with different directivities.
[27] Plastic deformation around a propagating shear crack can result in slower rupture acceleration and altered slip  accumulation compared with an elastically responding material with the same nucleation length and nucleation procedure, as shown in Figure 11. In an elastic material, the speed regime in which a shear crack propagates depends on the seismic S ratio [Andrews, 1976]. We have found that the transition to supershear in an elastic-plastic material depends on the extent and magnitude of plastic strain accumulation surrounding the fault, which are controlled by both the seismic S ratio (small values mean initial shear stress t xy 0 is already near peak strength) and CF (if CF ( 1, the plastic zone near the rupture front will be small).
[28] In Figure 11, for which the plastic response case corresponds to the large plastic zone in Figure 7a with CF = 0.93, S = 0.75, and Y = 14°, the supershear transition seems to be suppressed. However, in Figure 12, for CF = 0.61, S = 0.6, and Y = 56°, the supershear transition is only delayed rather than completely suppressed. For that case in Figure  12, the extent of off-fault plastic strain, shown along with contours of shear stress in Figure 13, increases with increasing rupture distance as the rupture accelerates toward the Rayleigh wave speed. After the rupture propagates a distance of 7R 0 , the shear wave peak traveling ahead of the crack tip reaches the peak strength of the fault and nucleates a daughter crack ahead of the main crack tip. Plastic strain accumulation begins around the daughter crack and the Mach front can be clearly seen in the shear stress contours. While the rupture propagates in the sub-Rayleigh regime, the maximum equivalent plastic shear strain, g eq p , close to the fault is 1.0t p /(2G), but the maximum equivalent plastic strain occurring along the part of the fault where rupture travels in the supershear regime is only 0.5t p /(2G). There is remarkable contraction of the plastic zone size during the daughter crack nucleation process.
[29] The propagation distance required for the transition to supershear depends on the seismic S ratio and the proximity of the initial stress state to elastic-plastic yield, CF. To gain some understanding of how CF controls the transition distance, we varied the fault friction parameters (with f d /f s = 0.1) so that CF could be controlled independently of S and Y. For Y = 56°, the propagation distance required for the transition to supershear at S = 0.75 increases from 9.6R 0 for CF = 0 (elastic off-fault material response) to 12R 0 for CF = 0.5 to 16.4R 0 for CF = 0.65 and to greater than 25R 0 for CF = 0.75. Rice et al. [2005] estimate R 0 to be on the order of tens of meters at midseismogenic depths for typically large crustal events, so 25R 0 will be a distance on the order of only a kilometer, and we cannot yet address what happens at greater lengths. Either S or CF must be Figure 15. Contours of equivalent plastic shear strain are shown for increasing levels of dilation. Nucleation lengths, L 0 c , and time after nucleation are identical for each case. For b ! 0.24m = 0.14 the critical hardening, h cr , as given by Rudnicki and Rice [1975], is nonpositive, and localization does not occur for h = 0. Here the grid dimension is Dx = R 0 /20. sufficiently low for the transition to supershear to occur within a finite propagation distance.

Plastic Strain Localization
[30] Bifurcation analyses for states of quasi-static deformation show that instabilities in the constitutive description for certain classes of elastic-plastic materials can lead to localization of deformation [Hill, 1952;Thomas, 1961;Rudnicki and Rice, 1975;Rice, 1976;Rice and Rudnicki, 1980]. The quasi-static localization conditions coincide dynamically with those for a vanishing propagation speed of elastic-plastic body waves [Hadamard, 1903;Hill, 1962;Mandel, 1963Mandel, , 1966Rice, 1976]. In the equivalent plastic strain fields of our dynamic analyses, tabular zones of high plastic strain whose cross section in the x À y plane has a needle-like appearance are notably visible for high angles of the most compressive stress, Y. In Figures 9a and 9b, the equivalent plastic strain fields for Y = 56°display parallel rows of the needle-shaped features. We show that these features are evidence of dynamic strain localization, occurring in the material due to pressure-dependent plastic yielding, and find that these elastic-plastic zones then fail to propagate perpendicular to their plane as waves. Instead, they grow in their own planes near the rupture front. The localization features become visible with increasing grid refinement, and their spacing is set by the element size. They are, as we show, largely understandable in terms of the existing theoretical background, and signal that no continuum solution actually exists for the adopted model. Some localization limiting procedure [Bazant et al., 1984;de Borst, 1988;Needleman, 1988;de Borst et al., 1999] would need to be added to the constitutive description in order for a solution to exist (for more details, see de Borst and van der Giessen [1998] and Bazant and Jirasek [2002]). That is an important goal for continuing work, although even without such a procedure we can gain some important perspectives. In particular, we now show that the conditions allowing the occurrence or not of localization are fully predictable from the existing theoretical background and that overall shapes of plastic zones and rupture propagation features in parameter ranges allowing localization are not very different from those for nearby parameter values which do not allow localization. While spacings at the smallest scale seem to be determined by grid size, our results suggest that the spacing of the longer localization features may in fact be independent of grid size and scale with a real physical length in the modeling, namely, with the slipweakening zone length R.  Figure 14). The grid size is Dx = R 0 /80.

Critical Hardening and the Presence of Strain Localization
[31] A pressure-dependent elastic-plastic material description, such as MC or DP, can be prone to localization even for positive hardening. Rudnicki and Rice [1975] found for the DP model that there is a critical hardening, h cr , above which localization cannot occur: Here m and b are the local slope of the yield surface and the plastic dilatancy factor, respectively, as defined previously, and N is a stress state parameter defined as where s 2 is the intermediate principal deviatoric stress. For this study, we use an initial stress state where the out of plane principal stress is equal to the average of the in-plane principal stresses. For that stress state, s 2 0 = 0 and N 0 = 0, so for a material with n = 0.25, h cr is positive (meaning a material with smaller h, like h = 0, would localize) when This condition is always met in our models neglecting plastic dilation (b = 0) and can be met when m > 0 for a range of sufficiently small (or large) b; b < 0.24m or b > 4.16m. (The latter range, with b > m, is unphysical in that it makes the plastic work rate s ij D ij p negative.) Equation (15) with N replaced by zero also describes the localization condition for the MC material, and equation (17) applies for that case too, at least in the nondilatant case, b = 0. If b 6 ¼ 0, we would have to consider if the present presentation of plastic dilatancy is what we would prefer in the MC case; if so, then equations (15) with N set to zero and (17) still apply.
[32] Locally, as preyield stress changes develop under plane strain, Figure 17. Contours of equivalent plastic shear strain for grid refinements of 20, 40, and 80 elements within R 0 , the static slip-weakening zone length, are shown for cases with zero hardening (h = 0) such that h < h cr 0 (= 0.031G). Nucleation lengths, L 0 c , and time after nucleation are identical for each case. The increasing prominence of localized plastic shear strain with increasing grid refinement shows that localizations may not become apparent except at extreme grid refinement, far beyond what is conventionally used in rupture dynamics studies.
So s 2 and hence N, become positive on the compressional side of the fault, stabilizing against localization in the DP material, but negative and hence promoting localization on the extensional side. Localization in the MC material is unaffected by changes in s 2 of equation (18), unless those changes are large enough that s zz ceases to be the intermediate principal stress. Figure 14 shows changes to the DP critical hardening ahead of the rupture tip for Y = 56°, with Dh cr = h cr À h cr 0 (where h cr 0 is the critical hardening based on the initial stress state) becoming positive and negative on the different sides of the fault, as anticipated. Changes in critical hardening promote plastic strain localization on the extensional side of the fault while stabilizing against localization on the compressional side.
[33] Including dilatant plastic deformation reduces h cr for a given m and n. For n = 0.25 and N = 0, as in the analyses presented above, the critical hardening h cr , will be nonpositive (no localization predicted) for b ! 0.24m. Figure 15 shows contours of equivalent plastic strain for Y = 56°with varying levels of dilation. When b ! 0.24m, h cr 0 0, and no localization features are present in the plastic strain field. The inclusion of substantial dilatancy (Rudnicki and Rice [1975] suggest b % 0.2 to 0.4 as representative for fissured rock masses) has a small effect on plastic zone size and propagation speed, at least in the present case of a ''dry'' granular material, without coupling to an infiltrating pore fluid as considered in part 2.
[34] When plastic strain hardening above the critical level is incorporated into the elastic-plastic description of the material, no strain localization is found in our numerical simulations, the result expected from theory. Figure 16 shows the equivalent plastic strains for Y = 56°with S = 1.0 and hardening ranging from 0.0h cr 0 to 2.0h cr 0 , based on the h cr 0 for N = 0. Localization features are completely absent when h = 2.0h cr 0 , which precludes even a local h cr > h, and that is accompanied by negligible alterations in the size and overall shape of the plastic zone, rupture propagation velocity, and total slip accumulation. We have therefore decided to leave for future work the far more elaborate gradient or nonlocal constitutive descriptions, and much more demanding computation challenges, which would allow resolution of localized structures. (Such would require grid sizes that are small compared to localization band thickness.)

Grid Dependence of Localization Features
[35] In the majority of our dynamic analyses we have specified no hardening (h = 0) and no plastic dilation (b = 0), so h < h cr 0 for the initial stress states considered with N 0 = 0. We observe shear-band-like structures of localized plastic strain at the scale of grid spacing, for grid spacings of R 0 /10 to R 0 /160. Figure 17 shows how the spacing of these localized strain features changes with changes in mesh resolution (Dx = R 0 /20 to R 0 /80) for Y = 56 with S = 1.0 and CF = 0.57. At high levels of grid refinement (Figure 18), the size and spacing of the longest localization features does not remain at the scale of grid spacing. Those features recur aperiodically along the rupture path, but have a characteristic spacing related, apparently, to the current length, R, of the slip-weakening zone, i.e, of roughly 0.2 to 0.4R, rather than to grid size. This leads us to the conjecture, to be tested in future, more computationally demanding simulations with a localization limiting procedure, that the spacing of the most prominent localization features scales with a real (i.e., grid-size-independent) physical length scale in the problem, namely, the slip-weakening zone size, R.
[36] Overall, the localization features are still strongly mesh-dependent. Nevertheless, for rupture propagation over distances on the order of 15 to 20R 0 like we have studied, the dynamics of rupture propagation (e.g., rupture length vs. time) does not vary significantly with mesh refinement, and seems to be little different from what is obtained when we increase the assumed plastic hardening modulus h above h cr Figure 18. A magnified view of contours of equivalent plastic shear strain is shown for Y = 56°for the same nonhardening (h = 0) case as in Figure 17 near the rupture tip for extreme grid refinements of (a) 160 and (b) 80 elements within R 0 , the static slip-weakening zone length, for rupture propagation distances of 9.5 R 0 and 12.6 R 0 , respectively. A physical (grid-independent) length scale, R, the current slip-weakening zone size, seems to scale the spacing between the longest localized strain features. Because of differences between Figures 18a and 18b in rupture propagation times, and hence lengths and speeds, the slip-weakening zone size is a different fraction of R 0 in Figures 18a and 18b. and obtain a locally smooth and presumably point-wise convergent numerical solution, as has been shown. The issue of localization must ultimately be confronted, both to provide a more rigorous basis for future studies of the type here, and especially because h must be expected to be negative in the early phases of deformation in a damage zone which has experienced partial cementation (hence gained a cohesion c) in the long period between earthquakes. That c would reduce toward zero with plastic strain in a new earthquake, meaning that h would be negative at small plastic strain, if perhaps transitioning to essentially zero as larger plastic strain develops. A sufficiently robust localization procedure would, of course, ultimately allow modeling of the spontaneous development of new shear fractures, and these new shear fractures should be expected to interact significantly with rupture propagation on the primary fault [Ando and Yamashita, 2007;Bhat et al., 2007b]. We have investigated localizations in material with zero or small hardening, however, in a softening material a wider spacing in the localization features may be expected because of a stress shadow effect of the localizations on the adjacent material.

Conclusion
[37] In this study we focus on the activation of distributed inelastic deformation during spontaneous dynamic shear rupture, in the crack-like mode, along faults that are bordered by zones of cohesionless damaged (cracked/granulated) material. We find that the location of inelastic deformation with respect to the fault depends strongly on the angle Y of most compressive stress to the fault of the initial stress field. In the present modeling for a ''dry'' porous material, inelastic deformation occurs only on the extensional side for higher angles of most compressive stress, Y ! 45°on both the extensional and compressional sides for 10° Y 45°, and only on the compressional side for Y < 10°. The extent of off-fault plastic deformation increases with increasing closeness CF of the initial stresses to the plastic failure threshold, and with decreasing seismic S ratio (of which small values indicate closeness of initial shear stress on the fault to failure). The presence of saturating ground fluids affects the location and extent of the zones of plastic activation as shown in part 2, because pore fluid pressure increases induced on the compressional side of the fault bring the material closer to violating the yield criterion, and suctions induced on the extensional side strengthen the material against yield.
[38] Off-fault plastic straining can have significant effects on the dynamics of rupture propagation and evolution of stresses along the fault. A change in the residual faultparallel stress, different on the two sides of the fault, occurs in the zone of plastic deformation, and this could be a signal for diagnosing rupture directivity. Further, the transition of rupture to the supershear regime in an elastic-plastic material is effected by the extent and magnitude of plastic strain accumulation surrounding the fault so that transition depends not only on the seismic S ratio, but also on how close the initial stresses in the material surrounding the fault are to yield levels.
Appendix A: Details of Elastic-Plastic Off-Fault Material Description [39] The Drucker-Prager yield function [Drucker and Prager, 1952;Lubliner, 1990, chapter 3.3.3] is given by where m and b may vary with plastic strain. A schematic of such a yield surface is shown in Figure A1. Plastic deformation first occurs when F(s ij ) = 0, that is, when t = Àm(s kk /3) + b. At least when h > 0, the yield function, F, can be used to determine whether an additional increment in stress, for a material at yield, will result in further plastic loading or elastic unloading: (@F/@s ij )ds ij < 0 for unloading, = 0 for neutral loading, and > 0 for loading. When h < 0, a test of sign of the numerator of the last fraction in equation (A7) to follow, which numerator is linear in the strain rates D kl , is necessary to distinguish plastic and elastic response; positive values (when the stress state is at yield) imply continuing plastic flow, a result that also applies when h > 0. Plastic strain increments can be obtained from the flow potential, M(s ij ) [see, e.g., Lubliner, 1990, chapter 3.2.1;Hill, 1950, chapter II.2 -4]. In the simple case that b is independent of the level of s, M(s ij ) is given by The plastic strain rate can be expressed as where Q ij (s ij ) and P ij (s ij ) are the gradients with resect to s ij of the yield function and flow rule, respectively: Strain hardening/softening occurs if b or m evolve with the plastic deformation, and the hardening satisfies hd p g/dt = _ b À _ ms kk /3 where ''equivalent'' plastic shear strain rate is (In the calculations which include hardening here, we fix m and let only b vary with plastic strain, _ b = hd p g/dt.) [40] Rudnicki and Rice [1975] note limitations of this type of description because simple modeling of cracked rocks suggests ''vertex'' effects so that P ij must have some dependence on the ''direction'' of _ ij . Also, a corotational measure of stress rate should more rigorously be used so that stress rate vanishes in rigid spin of material, but this effect is normally negligible for problems like those addressed here, where strains and rotations are small. This type of material description was shown by Rudnicki and Rice [1975], when b 6 ¼ m to permit localization of deformation for positive hardening, h, particularly in deformation states near plane strain.
[41] Total strain rates are given as the sum of elastic and plastic parts, as D ij = D ij e + D ij p , where D ij p is expressed above and We define L ijkl , with symmetries L ijkl = L jikl = L ijlk , such that equation (A6) inverts to _ s ij = L ijkl D e kl . In the present case of isotropic elastic response, L ijkl D e kl = (K À 2G/ 3)d ij D e kk + 2GD e ij , and that then fully defines L ijkl . The full elastic-plastic relations between stress rate and strain rate, assuming a stress state at yield and plastic loading, are Here the collection of terms after the minus sign corresponds to D kl p and the numerator and denominator in the last fraction of equation (A7) correspond, respectively, to Q pq L pqrs D rs and to h + Q mn L mnuv P uv . When the denominator is positive, as expected even for plausibly negative values of h, the condition for continuing plastic response is that the numerator be positive, as noted.

Appendix B: Split Node Contact Procedure
[42] We use a split node contact procedure to implement the friction law and to model unidirectional slip along the fault. Two contact surfaces define the fault, each with its own set of nodes, so that a duplicate set of nodes is defined along the fault, as shown in Figure B1, with nodal mismatch between the surface given by aDx, where 0 a 1/2. We use infinitesimal strain kinematics and consider changes in a due to slip to be negligible. The nodal momentum equations in the tangential (x) direction are given by where F j (1) and F j (2) are the nodal forces per unit distance in the z direction due to the stresses {s} in adjoining elements on surfaces (1) and (2).
[43] T j is the frictional force, per unit distance in the z direction, due to the surface interactions at each node, so that T j /Dx is the shear stress, t xy , along the contact surface. The nodal acceleration is given by x j and m j is the mass per unit distance in the z direction at each node along the contact surface. We use meshes with uniform node spacing so that m j is a constant, m. The slip velocity between the contact surfaces, linearly interpolated, has time rate À 2 T j þ að1 À aÞ T jþ1 þ T jÀ1 2 À T j ! ðB4Þ Figure A1. Yield surface (heavy line) for a material with pressure-dependent yielding with internal coefficient of friction, m and plastic dilation factor, b. Plastic strain increment components are superposed on the diagram. Since a (1 À a) < 1/4, and the quantity within the final parentheses should approach zero for good mesh refinement, _ V j can be approximated by where F j = F j (1) À (1 À a)F j (2) À aF jÀ1 (2) . For uniform time steps, integrating equation (B5) gives the slip velocity during the next increment: The force per unit distance in the z direction to prevent slip during the next time step is To implement the linear slip-weakening friction law (equation (2)) at each time step, the frictional strength of the fault, T j fric , is calculated at each node by where the T j norm are the nodal fault-normal forces, positive in compression, and f s , f d and D c are the static and dynamic coefficients of friction and the critical slip-weakening distance. This spilt note contact procedure is implemented within the ABAQUS/Explicit user subroutine, VFRIC. At each node, the minimum of T j fric and T j stop is applied as the frictional force per unit thickness between the contact surfaces, T j . Within the subroutine, T j stop is provided under the name f StickForce and the user must specify the frictional force between the contact surfaces, f Tangential, or T j .