Finite difference modelling of rupture propagation with strong velocity-weakening friction

SUMMARY

We incorporate rate- and state-dependent friction in explicit finite difference (FD) simulations of mode II dynamic ruptures in elastic media, using the Mimetic Operators Split-Node (MOSN) method, with adjustable order of spatial accuracy (second-, fourth- or mixed-order accurate), including an option that is fourth-order accurate at the fault discontinuity as well as in the elastic volume. At fault points, the rate and state equations combined with the spatially discretized momentum conservation equations form a coupled system of ordinary differential equations (ODEs) for slip velocity and state variable. As a consequence of the rapid damping of velocity perturbations due to the direct effect, this system exhibits numerical stiffness that is inversely proportional to velocity squared. Approximate solutions to this velocity-state system are achieved by two different implicit schemes: (i) a fourth-order Rosenbrock integration of the full system using multiple substeps and (ii) low order integrations (backward Euler and trapezoidal) of the velocity equation, time-staggered with analytic integration of the state equation under the approximation of constant slip velocity over the time step. In assessing the numerical schemes, we use three test problems: ruptures with frictional resistance controlled by (i) a slip evolution law with strong velocity-weakening behaviour at high slip rates, representing thermal weakening due to flash heating of microscopic asperity contacts, (ii) the classic (low-velocity) slip evolution law and (iii) the classic aging evolution law. A convergence analysis is carried out using reference solutions from a spectral boundary integral equation method (BIEM) (a method restricted to homogeneous media, with nominal spectral accuracy in space and second-order accuracy in time for smooth solutions). Errors are measured by root-mean-square differences of fault-plane time histories (slip, slip rate, traction and state). MOSN shows essentially the same convergence rates as BIEM: second-order convergence for slip and state-variable misfits, with slower (but at least first-order) convergence for slip rates and tractions. For a given grid spacing, fourth-order MOSN is as accurate as BIEM for all variables except slip-rate. MOSN-Rosenbrock nominally has fourth-order temporal accuracy for the fault-plane velocity-state ODE integration (compared to lower-order accuracy for the other two MOSN schemes) and therefore provides an important theoretical benchmark. However, it is sensitive to details of the elastic calculation scheme and occasionally its adaptive substepping performs poorly, leading to large excursions from the reference solution. In contrast, MOSN-trapezoidal is robust and reliable, much easier to implement than MOSN-Rosenbrock, and in all cases achieves precision as good as the latter without recourse to substepping. MOSN-Euler has the same advantages as MOSN-trapezoidal, except that its nominal first-order temporal accuracy ultimately leads to larger errors in slip and state variable compared with the higher-order MOSN schemes at sufficiently small grid spacings and time steps.

Rupture modelling with strong velocity-weakening friction 1833 in assessing (3). We address temporal discretization through numerical experiments with the time step size reduced below stability limits, following similar 3-D studies of  for the case of SW friction.
To sample a representative range of physical models, we vary the friction law and background stress load. We perform simulations with three different state evolution laws: the RS-FH law, as well as the classic Dieterich-Ruina aging (or slowness) law (RS-A) and (low-velocity) slip law (RS-S). We vary the background stress load so that both of the two well-known end-member modes of ruptures, expanding cracks and short-duration pulses, are included in the study. By varying these two features, state evolution and rupture mode, we are able to generate a range of test problems representative of recent theoretical, experimental and numerical studies (e.g. Zheng & Rice 1998;Lu et al. 2007).

In-plane ruptures under rate-and state-dependent (RS) friction
We model 2-D mode II dynamic ruptures on a frictional interface lying on the x axis. The fault line is surrounded by a linearly elastic and isotropic medium (x-z plane) with density ρ, Lamé constants λ and μ and P-and S-wave speeds c d = (λ + 2μ)/ρ and c s = √ μ/ρ, respectively. Frictional sliding is governed by RS constitutive laws expressed in the three different forms detailed below. The displacement and stress fields satisfy the plane-strain elastodynamic equations, with dots denoting time differentiation and commas denoting spatial differentiation, for example,ü x = ∂ 2 u x /∂t 2 , u x,x = ∂u x /∂ x, τ xx,x = ∂τ xx /∂ x, and so on. The components (u x, u z ) of the displacement vector and the incremental stress tensor (τ xx, τ zz and τ xz ) are measured with respect to initial distributions. The stress components τ zz and τ xz are continuous at z = 0, and for simplicity we also take u z to be continuous (although discontinuous u z , that is, fault opening, is easy to accommodate in the formulation). The fault-parallel displacement u x is discontinuous across the x axis and this jump is referred as slip s, that is, designating the limiting values for u x at this interface as u ± x (x, z = 0, t) = u x (x, z = ±ε, t) for ε → 0 + , then s = u + x − u − x . Friction laws used here are defined in terms of the magnitude of the time derivative of the slip, denoted by V = |u + x −u − x |, also referred to as sliding velocity. In the following, we use the superscripts (+) and (−) to denote limiting values at the plus (z > 0) and minus (z < 0) sides of the fault, respectively, of those wavefields and material properties that are discontinuous.
We adopt fault boundary conditions at z = 0 that allow spontaneous rupture propagation. For simplicity, we assume reflection symmetry of the elastic properties about the fault, with the consequence that the normal load, σ n , is time-invariant (the method generalizes easily to the asymmetric case with variable normal load). The total shear traction τ (x, t) acting on the minus side of the fault is composed of a known shear load τ load (x, t) and the incremental stress τ xz (i.e. τ = τ load + τ xz ). The magnitude of the total shear traction equals the frictional strength of the fault τ c (current value) and shares sign with the particle velocity discontinuity (but note V > 0, t ≥ 0), that is, The above represents a simplified case of more general formulations of shear rupturing (permitting, for example, stationary contact, normal stress fluctuations, fault opening) as in, for example, Day et al. (2005). Here, we assume that the frictional strength τ c is proportional to the normal stress σ n , with a proportionality coefficient f = τ c /σ n . The friction coefficient f depends explicitly on V and a variable θ that accounts for the state of the contact surface. We follow the usual (Dieterich 1978(Dieterich , 1979Ruina 1983) formulation for friction coefficient in terms of V and θ , which gives, Eq. (4) expresses the coefficient of friction as the sum of a reference value f * , the direct effect (giving the explicit velocity dependence, proportional to constant a), and the evolution effect (proportional to constant b). The aging law for state evolution iṡ where the state variable θ can be interpreted as the average lifetime of the asperity population and L as the average slip required for renewal of the contact population (see also, Beeler et al. 1994). In our numerical implementation, we favor the following equivalent expressions of (4) and (5) where The term f SS (V ) is the friction coefficient for steady-state sliding at fixed V, corresponding to θ SS (V ) = L/V , the steady-state value of θ in (5). The second RS evolution law we consider is the slip law, so called because state evolution occurs only during slip, with no evolution at V = 0 (Ruina 1983;Linker & Dieterich 1992). Transforming, as before, to dimensionless state variable ψ, the slip evolution equation iṡ with f SS given by (8). Note that (9) can be viewed as an approximation to (7), valid near the steady-state regime ( f ∼ f SS ). As a third evolution law, we use the slip law for state evolution (9), but with the steady-state coefficient (8) modified to account for strong weakening of friction at high sliding rates ( 0.1 m s −1 ). Models of f SS proposed by Rice (1999Rice ( , 2006 and Beeler and co-workers (Beeler & Tullis 2003;Beeler et al. 2008) to simulate flash-heating of asperities include a critical velocity V w , above which f SS decays rapidly, while retaining the weak (logarithmic) velocity dependence observed in laboratory experiments at lower slip rates. Specifically, we combine (9) with a steady-state friction coefficient given by Estimates of f w (fully weakened coefficient of friction) are in the range 0.2-0.3 (e.g. Prakash 1998Prakash , 2004; and references therein), and estimates of V w (the weakening velocity) are in the range 0.1-0.5 m s −1 (Rice 1999(Rice , 2006. We abbreviate the three cases as RS-A (for the aging law, eqs 6-8), RS-S (for the slip law, eqs 6,8,9) and RS-FH (for the enhanced weakening law, eqs 6, 9, 10).

Test cases: parametrization, rupture nucleation and initial conditions
We aim to quantify accuracy of our numerical algorithms under the three constitutive laws, while illustrating features of ruptures that emerge under these friction laws at high rupture velocities and high slip rates (V ∼ 1 m s −1 ). We generate numerical test cases for each evolution law, using a common reference set of homogeneous medium and frictional parameters (described in Table 1), with a few exceptions that we will clarify. In particular, the slip length scale L varies according to the friction law: L = 0.06 m in the case of RS-A, while L = 0.4 m in both RS-S and RS-FH. Our use of different L values was motivated by the expectation that the effective slip-weakening displacement for rupture simulations with the aging evolution law would be larger than for those with the slip law when the same L is used with both evolution laws (e.g. Bizzarri & Cocco 2003). We used a larger L for slip law simulations in order to achieve comparable effective slip-weakening displacements for all test problems (although, as discussed later, our simulations still produced somewhat larger effective weakening displacements for RS-A compared with RS-S and RS-FH). The shear load, τ load (x, t), acting on the fault at any time is the sum of a constant background stress, τ 0 , and a perturbation term, τ p (x, t), that is, τ load = τ 0 + τ p . The perturbation term is used to nucleate 2-D mode II ruptures, and takes the form, This mode of nucleation is not physically motivated, but rather is intentionally designed to nucleate a smooth solution, in order to facilitate study of convergence rates of the numerical schemes employed. Note that τ p is Gaussian in space, with standard deviation w, and increases smoothly over time interval [0, T ] from zero to its final amplitude A 0 . All derivatives of (12) are continuous and of finite support, making this function useful for testing high-order numerical schemes (i.e. this loading function does not introduce discontinuities, of any order, that might mask the nominal convergence rate of a scheme). Noda et al. (2006Noda et al. ( , 2009) used a model similar to eq. (11), but step-like in time, to nucleate mode III ruptures. Table 1 lists the values for σ n , τ 0 , A 0 , w and T that are used in our simulations, and shows that choices for (τ 0 , A 0 ) depend on the evolution law. These τ 0 values result in crack-like solutions in test cases for RS-A and RS-S, and pulse-like ruptures for RS-FH, permitting the numerical solvers to be tested for both rupture modes. Finally, we assign A 0 such that peak slip rates are typical of coseismic values. We assume uniform initial conditions in displacements u, velocitiesu, and the state variable ψ, with u = 0 andu z = 0 everywhere for t ≤ 0. Along the fault, the initial sliding velocity V 0 > 0 is uniform, anḋ Condition (13) makes the fault line a symmetry (or antisymmetry) axis for elastic wavefields. Initial state ψ 0 must satisfy the boundary condition (3) subject to an initial sliding rate V 0 , and we accomplish this by using the integral form of the elastodynamic equation along z = 0 (Cochard & Madariaga 1994;Cochard & Rice 1997), In (14), τ load (x, t) is the shear load (i.e. stress on the fault in the absence of slip), (2c S ) −1 μV is the radiation-damping response (which relates instantaneous changes in stress to those in slip rate), and φ(x, t) is the stress transfer functional (Cochard & Madariaga 1994;Fukuyama & Madariaga 1998), which depends upon fault slip, vanishing when slip is uniform. Since τ load = τ 0 + τ p , and both τ p (x, 0) = 0, and φ(x, 0) are zero at t = 0 eqs (3) and (6), leads to Eq. (15) enforces an initial state ψ 0 consistent with V 0 and τ 0 . We take V 0 = 10 −12 m s −1 (∼0.03 mm yr −1 ) as fixed in our three test cases of constitutive equations.

Solution assessment
Our objective is to test an explicit FD method, based on the MOSN algorithm of Rojas et al. (2008), for modelling rupture under each of the three RS laws. MOSN uses a split-node technique to incorporate the fault discontinuities, and Rojas et al. (2008) previously applied it to the case of purely slip-dependent friction. We assess the new method by comparing solutions with reference solutions calculated with a spectral BIEM of Noda et al. (2006Noda et al. ( , 2009. In these comparisons, we adopt the same convention introduced in Day et al. (2005) to distinguish individual simulations according to both the method and the grid size h (measured in metres) used. Thus, calculations for h = 40 m, for example, are denoted by BIEM040 and MOSN040, respectively, for the BIEM and MOSN methods. To be consistent with the 1-D discretization of the fault line performed by BIEM, we only compute MOSN simulations on uniform square grids, that is, a common grid size h is used in both spatial directions, x and z. Given that no analytical solutions are known for the test problems, we use as our reference solution in each case a relatively fine-meshed (h = 10 m) BIEM solution, BIEM010. Reference solutions BIEM010 will be used in: (i) evaluating alternative implementations of MOSN boundary treatment at fault points, (ii) estimating the transient size of the cohesive zone observed in test problems under consideration and (iii) assessing the convergence of MOSN and BIEM solutions under grid refinement.

M I M E T I C O P E R AT O R S A N D S P L I T N O D E S ( M O S N ) M E T H O D
The MOSN FD method is described in detail by Rojas et al. (2008). Medium properties and elastic fields are discretized using a rectangular 2-D staggered grid augmented with split nodes placed at the fault plane (see Fig. 1). Apart from this boundary, this grid reduces to that of Madariaga (1976) and Virieux & Madariaga (1982), and is sometimes referred to as a Standard Staggered Grid (SSG) (e.g. Bohlen & Saenger 2006). Following the method of , the split nodes represent both stress and displacement (and velocity) discontinuities on the fault, in the process introducing certain auxiliary stress and displacement discretization points. The main innovation of MOSN was the introduction of a consistent fourth-order accurate FD discretization of these fault-plane discontinuities, achieved by using mimetic FD operators. Optionally, artificial viscosity of Kelvin-Voigt form can be used to attenuate any spurious high-frequency oscillations (excited by, e.g. the discrete advance of the rupture through the numerical lattice). Split nodes also accommodate extra evaluations of shear stress τ * xz and normal displacement u * z to permit the approximation of τ xz,z and u z,z using the mimetic operator G. These approximations, in addition to values of slip rate V and state variable ψ at u x − sites, are used by our TSN implementation of faulting under rate-and state-dependent friction (eqs 7-10 and 19-20). Grid locations of material properties (ρ, μ and λ) correspond to original definitions given by Madariaga (1976), while interpolations of densityρ and shear modulusμ follow interior averaging strategy of Graves (1996) and fault-plane approximations of .
For the current study, we add two main features to MOSN: (i) rate-and-state friction with the three evolution laws (eqs 7-10), including implicit ODE solvers for slip velocity and state variable (with a posteriori computation of traction and fault-slip) (as one case of which was discussed by Rojas et al. 2007) and (ii) options for reducing the order of accuracy of the spatial differentiation operators (to facilitate a more complete analysis of numerical convergence and accuracy). The resulting numerical tests will thus also be relevant to the potential performance of RS implementations in other (2-D or 3-D) explicit split-node FD, finite element, and support operator codes, most of which are either second order (e.g. Oglesby et al. 1998;Aagaard et al. 2001;Bizzarri & Cocco 2005;Day et al. 2005;Ma & Archuleta 2006;Ely et al. 2009) or mixed order (e.g. , Moczo et al. 2007a).

Second-, fourth-and mixed-order implementations
Here, we use FD approximations of second-, fourth-and mixed-order accuracy to approximate spatial derivatives of the system (1)-(2). The second-and fourth-order FD formulas each maintain consistent order of accuracy everywhere in the domain (including the fault plane discontinuity), while the mixed-order case is fourth-order accurate except in the approximation of z-derivatives at fault gridpoints (where it reduces to second order). The fault-plane discontinuity of certain wavefields (u x , u z , τ xx ) is accounted for in approximating z derivatives by using the mimetic differentiation matrices D and G given in Appendix A, which reduce to one-sided stencils at the fault plane (where necessary, subscripts are used to distinguish the three order-of-accuracy cases, e.g. D 2,4 incorporates second-order fault-plane stencils and fourth-order interior stencils). Fig. 1 shows the positive side (z ≥ 0) of the 2-D SSG, as in Rojas et al. (2008), but with the addition of slip rate V and state variable ψ on fault nodes (doubled line in Fig. 1 coincident with the x axis). Locations of gridpoints are given by (x, z) = (i, j) h where each of these indices takes values in the set {0, ±1/2, ±1, ±3/2, ±2, . . .}. We assign integer pairs (i, j) to grid sites of displacement components u x , and indexing of other discrete quantities (wavefields and material properties) is given by location relative to (i, j) in the unit cell. At fault nodes (indexed by j = 0), both the material properties (ρ, λ, μ) and the discontinuous field variables (u x , τ xx ), are assigned separately to each parts of the split nodes, and denoted as ρ ± i , λ ± i+1/2 , μ ± i+1/2 , (u ± x ) i and (τ ± xx ) i+1/2 (omitting the j index in the presence of plus or minus superscripts). Similar notation is used for the discrete fault variables V i and ψ i , and for auxiliary fault-node evaluations of vertical displacement (u * z ) i+1/2 and stress component (τ * xz ) i (see Rojas et al. 2008, for definitions of the latter quantities). The MOSN methodology can accommodate any scheme for the assignment of material properties to the indexed points, such as the volume averaging method of Moczo et al. (2002). Fig. 1 illustrates, as an example, the approach of Graves (1996), in which material properties are assigned at one point per unit cell, and Rupture modelling with strong velocity-weakening friction 1837 then interpolated as needed at other points (with interpolated density,ρ, obtained by arithmetic averaging, and interpolated shear modulus μ, obtained by harmonic averaging). The scheme is modified, as in , to give one-sided properties averages at the split nodes [see eqs (18) and (19) in Rojas et al. 2008].
The numerical procedure for z differentiation is more laborious and follows Rojas et al. (2008). It requires a special arrangement of discrete fields u x , u z , τ zz , and τ xz (including the fault variables (u ± x ) i ,(u * z ) i+1/2 ,(τ zz ) i+1/2,0 and (τ * xz ) i ) into two separate vectors v and f, in one of two forms, depending upon the x index. In the case of gridlines indexed by constant (integer) values i of the x index (containing the split u ± x sites), vectors v and f are defined as Approximations to u x z and τ xz z at required sites along this gridline are computed by products, Dv and Gf, respectively (see Appendix A, and note that the current definitions of D and G differ by a factor 1/h from those of Rojas et al. 2008). At gridlines indexed by constant (half-integer) values i + 1/2 of the x index (containing the split τ ± xx sites), vectors v and f are defined as and τ zz z and u z z are approximated by Dv and Gf, respectively. The displacement field at interior gridpoints ( j = 0) is time stepped by an explicit scheme based on centered, second-order discretization of the time derivatives in eq. (1), and the scheme is stable for time step satisfying t 0.6h/c d . At the split nodes, the equations of motion (1) are coupled to the friction eqs (6)-(10), and discretization and time integration of those coupled equations are discussed in the subsequent subsections.

Traction at Split-Node (TSN) equations: velocity-state ODE system
We use the approach of Rojas et al. (2008) to couple frictional boundary conditions (3) to the equation of motion (1.1), with modifications to accommodate RS laws (eqs 6-10). An extensive discussion and justification of the following numerical method is given in Appendix B in the simplified context of 1-D elasticity. The semi-discrete (i.e. spatially discretized) version of eq. (1.1), in terms of velocities v ± x at split node (i, 0), is where g 1j are components of matrix G K (a subblock of differentiation matrix G) given in Appendix A (eqs A6-A8). The stencil (g 11 , . . . , g 16 ) corresponds to a forward stencil appropriate to approximate τ xz z on the plus-side of (i, 0), while the minus-side approximation to τ xz z is obtained by using the backward stencil (−g 16 , . . . , −g 11 ). Specific entries g 1j depend on the order of accuracy of the simulation under consideration. Following Rojas et al. (2008), we rewrite (18) in the form of a force balance. First, we identify split nodal masses M ± i = ρ ± i h 3 /2, and the restoring forces R ± i (acting on plus and minus sides, respectively) that correspond to interior stress terms in (18) (i.e. excluding the fault shear-stress term τ * xz ). Next, the term τ * xz is written as τ * xz = τ − τ load by using the definition of total shear traction τ given earlier. Then, using the boundary condition (3) to express τ in terms of fault strength τ c , eq. (18) takes the form, where Eq. (19) reflects the dependence of split particle accelerationsv ± x on slip rate V (t) and the state variable ψ(t) through the constitutive law (6) for fault strength τ c . This dependency of τ c on ψ couples (19) to evolution eqs (7) and (8) for RS-A, (8) and (9) for RS-S or (9) and (10) for RS-FH. Hereafter, we call these coupled ODEs the velocity-state system and explain below two alternative integration approaches. The first is a relatively high-order scheme and thus provides a theoretical benchmark that is valuable in helping isolate effects of the spatial discretization. The second (which itself has two variants) is nominally of lower order of accuracy, but is much simpler to implement, entails no loss of accuracy in practice and proves much more reliably stable than the first method. x and state variable ψ are coincident and defined half t from stresses (τ ± xx , τ zz , τ * xz ), and restoring forces R ± , the distribution used in the Rosenbrock integration of the coupled velocity-state system (eqs 7-10 and 19-20). In panel (b), split velocities v ± x are half-t staggered relative to state variable ψ and stresses (τ ± xx , τ zz , τ * xz ), showing the integration of the velocity eq. (19) through t by backward Euler or trapezoidal rules. In both cases, shear traction τ is obtained from velocity through boundary condition eq. (6), split tangential displacements u ± x and normal displacements u * z are defined simultaneously with all fault-plane stress components, and fault-plane differential stress τ * xz is obtained from total traction τ (τ = τ load + τ xz ).

High-order unstaggered integration of the velocity-state system
Figs 2(a) and (b) compare two alternative time distributions of the fault variables. Both time-stagger velocities v ± x (by t/2) with respect to coincident displacements, stresses, traction τ and restoring forces R ± i . The fact that stresses are simultaneous to displacements and one-half time step staggered from velocities is consistent with standard velocity-stress formulations (e.g. Virieux 1986). Figs 2(a) and (b) differ in the time location of the state variable with respect to velocities. In Fig. 2(a), variables ψ and v ± x coincide and are time-stepped simultaneously by a high-order stiff ODE solver (detailed below). In Fig. 2(b), in contrast, variables ψ and v ± x are time-staggered and in this case we integrate slip-velocity eq. (19) by means of a low-order implicit scheme, while the state variable is updated semi-analytically, as explained in the next subsection.
The high-order time integration scheme splits into three steps.
Step 1 calculates both displacements, i are known from previous calculations (or from initial conditions), a centred integration yields u ± Next, we calculate [u * z ] into Hooke's law (eq. 2), differenced to get the split stress gradient t 0 + t i by integrating simultaneously the state-velocity system (see Fig. 2a) by numerical means. This integration requires a routine designed for stiff ODE systems in order to successfully handle frictional parameters and initial conditions considered in our simulations (see, Appendix B). We use a Rosenbrock algorithm of order (3)4, as given by Hairer & Wanner (1996) (code RODAS, available at www.unige.ch/hairer/software.html). This method is a fourth order scheme that includes an embedded third order scheme; the difference in the values obtained by the two schemes provides an estimate of the numerical error. The error estimate is then used to adapt the step sizes to bound error below a specified tolerance. Integration on [t 0 , t 0 + t] yields half-way values of velocities [v ± x ] t 0 + t/2 i (from which we form slip rate V t 0 + t/2 i ) and state ψ t 0 + t/2 i , then total traction τ t 0 + t/2 i is obtained from constitutive law (6) followed by boundary condition (3), and finally τ * xz = τ − τ load . Rosenbrock integration substeps each split node through the time interval [t 0 , t 0 + t], requiring the evaluation of right-hand terms of (19) at instants t, t 0 ≤ t ≤ t 0 + t. The restoring forces R ± , however, depend upon off-fault displacements. One of our principal goals is a scheme that will work with existing explicit elastodynamic codes, so we avoid coupling the fault-plane time integration to the time integration of off-fault field variables. With this constraint, R ± values are only available at every half-integer time step t 0 + i t/2 (i = . . . , − 3, − 1, 1), so we require some interpolation procedure. Rojas (2009) tested some alternative interpolation strategies in previously detailed RS-FH test (see Table 1) from which the following linear approximation gives adequate results, In Appendix B, we quantify the stiffness S of the above state-velocity system by examining the eigenvalues of the linearized system, finding an inverse dependence of S on V 2 . By stiffness, we mean that the set of ODEs contains variables that evolve over vastly different timescales, and we define S as the ratio of the two timescales in the velocity-state system. Explicit integration methods must resolve evolution of the fastest variables for stability, which often vary over scales that are much smaller than those of interest. Thus, explicit integrators become incredibly inefficient since they must take a possibly huge number time steps to span the much longer time scale of interest. An implication of this result is that under initial slip rates V 0, typical of interseismic regimes (e.g. V 0 ∼ 10 −12 m s −1 ), this system becomes highly stiff (S ∼ 10 −24 ) and even various standard stiff integrators fail in early stages of such simulations. In particular, we compare the Rosenbrock method with two alternatives tested: MEBDFDAE (based on backward differentiation formulae; Cash 2000), and RADAU5 (a fifth-order Radau Runge-Kutta method; Hairer & Wanner 1996). MEBDFDAE and RADAU5 succeeded only after initial conditions are restricted to sliding velocity V 0 > 10 −12 m s −1 . Bizzarri et al. (2001) and Bizzarri & Cocco (2005) also pointed out the utility of Rosenbrock's technique in FD rate-and-state calculations (without enhanced velocity weakening). Nonetheless, the Rosenbrock integration of the velocity-state system is not very robust, in the sense that it is quite sensitive to details of the elastic calculation scheme and the time and spatial step sizes. In some cases its adaptive substepping performs poorly, leading to large excursions of the Rosenbrock solution from the reference solution (but without clear indicators of instability). We have found solution anomalies of this type to occur due to small changes in the restoring-force interpolation scheme, changes in the elastic grid from staggered to partly staggered (as defined by, e.g. Moczo et al. 2007b), and small increases in grid spacing or time step. The alternative scheme introduced in the next subsection does not present these difficulties for any test considered here.

Staggered integration of the velocity-state system
This approach retains the same steps 1 and 2 as before, but the velocity-state integration, step 3, differs. With velocity and state time-staggered as in Fig. 2(b), we integrate the state evolution eqs (6)-(10) over the interval [t 0 − t/2, t 0 + t/2], with slip velocity V approximated by its value at the center of that interval, t 0 . With this approximation, the state equations become linear ODEs that we can integrate analytically, and in the particular case of slip evolution (eqs 7-10) the result (at fault node with index i) is given by An analogous integration can be applied in the case of RS-A evolution modeled by eqs (6) and (7).
In the final step, we integrate the slip velocity V from eq. (19) over the interval [t 0 , t 0 + t], now approximating ψ by ψ t 0 + t/2 . We have found that either of two low-order implicit schemes, backward Euler or trapezoidal, works very reliably, yielding a nonlinear algebraic equation that is easily solved numerically. For the backward Euler case, for example, we discretize eq. (19) with known termsṽ ± i given bỹ Note that eqs (24) and (25) involve evaluations of the state variable ψ and restoring forces R ± available at t 0 + t/2 which are time-staggered with respect to the discretization point; hence, no time interpolation is necessary. Next, we use eq. (25) to write the slip rate component v A crucial point of this formulation is that slip rate component v shares sign with ṽ i under the physical constraint of positive coefficient of friction (see eq. 6) and the fact that c < 0 (g 11 is also negative, see Appendix A). Then, eq. (26) can be written as an algebraic equation in [V ] t 0 + t i after multiplying by sign( v t 0 + t i ). This new equation is non-linear due to the constitutive definition of fault strength (6) whereṼ i = | ṽ i |. Newton's method solves eq. (28) quickly, especially after (28) is rewritten in terms of dimensionless variable w = ln(V i /V * ) serves well as the initial guess. Upon convergence, we invert this change of variable and the solution [V ] t 0 + t i is plugged into eq. (24) to yield individual split velocitiesṽ ± i . A trapezoidal rule version is nearly as simple to implement, is nominally of higher (second-order) accuracy compared with backward Euler, and gives superior global precision of the MOSN solutions, as shown in next sections. The trapezoidal integration introduces minor changes of above formulation and reduces to a substitution of the factor t by t/2 in eqs (24), (26), and (28), and a redefinition of the auxiliary variableṽ ± i by the following equation (instead of 25) The three alternative MOSN integration schemes for the velocity-state system given in this section, high-order Rosenbrock, time-staggered Euler, and time-staggered trapezoidal, can each be modified to account for 2-D sliding problems with variable normal stress σ n (not considered in current tests). In such cases, integration (step 3) must be preceded by the interpolation of [τ zz ] i+3/2,0 ], in cases of second-order simulations. An alternative fourth-order averaging formula for [τ zz ] t 0 + t/2 i,0 would be appropriate for higher-order calculations. In the rest of this paper, we refer to these three different MOSN schemes as the abbreviation 'MOSN' preceding the ODE method of numerical integration at fault points. That is, 'MOSN-Rosenbrock' names the unstaggered high-order scheme described above, and similar naming is used in cases of Euler and trapezoidal staggered integrations.

C O H E S I V E Z O N E E S T I M AT I O N
In this section, we estimate the size of the cohesive zone in each of our test problems. The cohesive zone is the fault region where shear traction drops and slip accelerates in response to the dynamic weakening process taking place behind the crack tip, and the extent of this high-gradient zone provides an estimate of the minimum length scale to be resolved by the numerical simulations. To measure , we study the behaviour of fault variables (traction, sliding velocity and state) in a neighbourhood of the crack tip and identify aspects of the breakdown process that lead to a practical definition of the cohesive zone applicable to the test problems under consideration.
To illustrate the evolution of fault variables within the cohesive zone in a single plot, we express these variables in dimensionless terms. We divide the integral form of the elastodynamics eq. (14) by the normal stress on the fault σ n (a constant parameter in this work), that is, The dimensionless fields of shear traction and slip rate are given, respectively, by the left-hand term and the second right-hand term (omitting the minus sign) of (30). We denote ( We therefore introduce an alternative cohesive-zone estimate. At a given time t = t 0 , we construct the tangent line with maximum slope to the curveτ (x, t 0 ) on x f < x < x u . We then define x p as its intersection with the constant levelτ f (see Figs 3a-c) and use the interval [x p , x u ] as a practical estimate of the cohesive zone. Figs 3(a)-(c) confirm that, in all three cases of the evolution law, most of the degradation of shear traction, and most of the state evolution, occur in the interval [x p , x u ], as does the maximum slip rate. In addition, this interval captures most (∼75 per cent on average) of the slip rate increase at the rupture front (i.e. the slip rate at x u is only about 25 per cent of its peak value, based on averages of 500 BIEM010 records equally spread on the time interval 3 ≤ t ≤ 8 s). We adopt (x u − x p ) as our estimate. Fig. 4(a) shows the locus of points x f , x p and x u for time 0 < t < 12 s for the three evolution laws (from reference solution BIEM010), with cohesive zone size corresponding to (x u − x p ), and Table 2 lists the minimum ( min ) and median (¯ ) size in the fault interval [4,15] km (avoiding the nucleation zone). We use these values to express numerical errors in terms of a dimensionless measure of numerical resolution, h/¯ (or, h/ min ), in the next section. Note the large effect of the state evolution model on , with case RS-A yielding magnitudes of min and¯ nearly 3.5 times bigger than those of RS-S, and more than 4.3 times bigger than estimates from RS-FH.    Fig. 4(b) shows the implied slip-weakening curves at the fault location x = 6 km. Note that curves for RS-S and RS-FH are very similar, and have slightly higher peak tractionτ u ; larger, more non-linear weakening rates; and smaller implied weakening distances d 0 , compared with the almost linear weakening curve for RS-A. We apply the same tangent-intercept method as before to estimate d 0 from SW curves, our d 0 being the slip value at which the tangent line with maximum weakening rate equals the minimum traction. Table 2 lists the median valuē d 0 (on fault interval [4, 12] km for BIEM010) for each test case. Values are nearly identical for RS-S and RS-FH, but the result for RS-A is about a factor of three larger.

E R RO R A S S E S S M E N T O F N U M E R I C A L S O L U T I O N S
Using BIEM010 as reference, we examine the convergence rates of MOSN for our three test cases, compare with those of BIEM, and quantify the effects of: (i) artificial damping (with parameter η), (ii) order-of-accuracy of spatial differentiation and (iii) time discretization. We express grid resolution in terms of the number of gridpoints discretizing the length scale¯ for a given h, that is,N c =¯ / h. For each numerical  Table 3 lists the corresponding values of h,N c and the time step t (setting t = 0.5 h c −1 d , giving Courant-Friedrich-Lewy number cfl = 0.5).

Global metrics
As our estimate of the numerical error in time-dependent quantity s(t) from a test simulation, we define a quantity MAX-RMS. MAX-RMS gives the misfit of time-series s i relative to the corresponding series s REF where the sum in (31) is over all discrete time and the maximum is taken over all fault points in the segment [0, 15] km. Metric (31) serves as a global measure of accuracy, accounting for misfits in both space and time. The reference grid (h = 10 m) contains all other grids used in the simulations, which avoids any requirement for spatial interpolation in applying metric (31).

Artificial viscosity
With an imperfectly resolved cohesive zone, the discrete advance of the rupture front through the grid may excite artificial high-frequency waves which contaminate the solution. Kelvin-Voigt viscosity has been successfully used in volume-discretized rupture simulations to mitigate this numerical artefact (e.g. Day 1982Day , 2005Oglesby et al. 1998Ben Jemaa et al. 2007). Rojas et al. (2008) used artificial viscosity to reduce numerical artefacts in simulations with MOSN using slip-weakening friction, and we here explore this approach to RS simulations. Elastic stresses (eq. 2) are modified by stress-rate-proportional damping terms to give modified stresses [τ xx ] i+1/2, j , [τ zz ] i+1/2, j , and [τ xz ] i, j+1/2 of the form (componentsτ zz andτ xz have similar definitions), andτ replaces τ in the equations of motion (and therefore also in the split restoring force terms R ± ). In (32), the coefficient η (η > 0) is a dimensionless damping parameter. Given that t is proportional to the grid size h ( t = 0.5 h c −1 d ), the absorption-frequency spectrum for constant η scales with h, and only wavelengths near the grid Nyquist limit are significantly damped (Day et al. 2005). We have confirmed through numerical tests that MOSN, with each of the proposed velocity-state integration schemes, performs well for all relevant values of η.
We examine solution sensitivity to η in problem RS-FH (which is representative of results from the other two test problems) only in the case of the MOSN-Rosenbrock method. Fig. 5 shows curves of MAX-RMS for slip rate and state as functions of η and h, with fourth-and second-order slip-rate solutions in parts (a) and (b), respectively, and fourth-and second-order state solutions in parts (c) and (d), respectively. For fourth-order MOSN, non-zero η provides only negligible improvement in slip rate (5a), at the cost of significant degradation in the state variable calculation (5c), so we conclude that η = 0 is the preferred value for the fourth-order algorithm. For second-order MOSN, η = 0.025 usually provides some very minor improvement in slip rate, and sometimes in both variables (5b, d), but this comes at the cost of equally minor degradation of the solution for traction and integrated slip (not shown). So for the second-order case, η values up to 0.025 are admissible, but there is no clear advantage to the use of non-zero η (we examine this issue further in the next section). We also note that, in every case, the dependence of MAX-RMS on η diminishes with decreasing h, becoming, as expected, negligible as the cohesive zone becomes better resolved (as also noted by Day et al. 2005).
The conclusion that artificial damping has little or no advantage for RS problems with MOSN differs from our finding for MOSN simulations of rupture under slip-weakening friction (Rojas et al. 2008). There we found significant improvements with the use of artificial viscosity, and other authors have found similar improvements when using artificial viscosity with slip weakening models (e.g. Day et al. 2005;Figure 5. MOSN misfits (MAX-RMS metric, eq. 24) as function of the artificial viscosity η for different grid sizes h. The test problem corresponds to RS-FH, with solution BIEM010 as reference. Misfits in slip rate and state variable, respectively, are shown in panels (a) and (c) for the case of fourth-order simulations, and in panels (b) and (d) for the second-order case. Dalguer & Day 2007, for FD methods;Ben Jemaa et al. 2007, for a finite volume method). We attribute the better performance of undamped RS models, compared with undamped slip-weakening models, to the explicit dependence of fault strength on slip velocity (i.e. the direct effect) that is built into the RS models (in agreement with experimental data). The explicit velocity dependence acts as a velocity-strengthening term that opposes slip acceleration (since the constant a in eq. 6 is always positive), and velocity weakening enters only implicitly, through the state evolution.

Convergence under grid refinement
Time-differencing accuracy in MOSN is nominally second order (except in the case of MOSN-Euler, which formally reduces to first-order due the first-order accuracy of the on-fault velocity-state integration). Thus, when the time step t is tied to the spatial step h through the CFL constraint t = c f l h c −1 d (as is conventionally the case), there is the potential (in spatially fourth-and mixed-order MOSN) for high-order convergence in h to be masked by lower-order errors due to time discretization. The latter effect, if it were present, might call for auxiliary schemes to equalize space-and time-discretization errors, for example, Lax-Wendroff corrections (for which seismic modelling examples are given by Blanch &Robertson 1997 andChen 2007) or integration with a high-order explicit Runge-Kutta method. Following , we perform some additional tests to assess whether apparent convergence rates of fourth-order MOSN (with respect to h) are contaminated by time-differencing errors by artificially reducing the CFL number below its stability limit. That is, we perform the coarsest-gridded MOSN simulation with steps (h 0 , t 0 ) satisfying CFL number cfl 0 = 0.5, and perform finer-gridded simulations with CFL In the latter case, calculations oversample the time axis by reducing the Courant-Friedrich-Lewy number cfl below its stability limit (∼0.6) as h decreases (following eq. 33). All other cases use a constant CFL number clf = 0.5. number c * (h) proportional to h, that is, As a result (and for these additional tests only), the time step is proportional to h 2 ( t ∝ h 2 ), and time-differencing errors scale with h 4 . We first use problem RS-FH and the MOSN-Rosenbrock method to assess convergence rates, under the MAX-RMS metric (relative to BIEM010), for all on-fault variables (sliding velocity V , traction τ , slip s and state-variable ψ). The analysis is done for each of the three spatial differencing schemes, and the contribution of artificial viscosity (eq. 32) is explored in second-order calculations. In addition, we evaluate the effect of low-order time differentiation by means of the above device (eq. 33). Fig. 6 shows MOSN-Rosenbrock errors as a function of grid size for sliding velocity V (Fig. 6a), traction τ (Fig. 6b), slip s (Fig. 6c) and state-variable ψ (Fig. 6d). The abscissa is the median cohesive-zone resolutionN c, and labelling along the top shows the corresponding grid size h (see also Table 3). For comparison, BIEM results are also shown. MOSN solutions follow approximate power laws, with errors in state variable (Fig. 6d) as the best expression of this tendency. BIEM results also follow power laws if we discount outlier BIEM060. We estimate the order of accuracy by linear least-square fitting on the logarithmic scale, and results are given in Table 4 (we include BIEM060 in this estimation). For reference, the solid line in each panel of Fig. 6 shows quadratic (second-order) convergence.  Fig. 6 shows very similar accuracy for fourth-and mixed-order MOSN. The fourth-order simulations are slightly more accurate in modelling shear-traction (e.g. ∼34 per cent more accurate at h = 50 m), but misfits from the two methods are nearly indistinguishable for slip-rate, fault-slip and state-variable misfits, especially for grids with h ≤ 50 m (N c ≥ 5). These schemes display approximately linear convergence in terms of slip rate and traction, and approximately quadratic convergence for fault slip and state variable (Table 4).
Second-order calculations are in all cases less accurate than higher-order MOSN solutions in low-resolved grids (h ≥ 50 m,N c ≤ 5). In that range of h, misfits in slip rate, fault slip, and state variable from second-order calculations are at least 100 per cent higher than the corresponding higher-order results (up to 300 per cent for state-variable errors). Similar error excess is just 50 per cent (also, h ≥ 50 m) on traction misfits. For the state variable, this accuracy gain (for fourth-and mixed-order over second-order) persists also for the better-resolved grids, since all three methods share a roughly quadratic convergence rate. For the other fault variables, however, second-order MOSN shows somewhat higher convergence rates than the nominally higher-order schemes, and overcomes much of the advantage of the latter for the most well-refined meshes. Fig. 6 also demonstrates that second-order time discretization is not a significant factor for MOSN-Rosenbrock accuracy or convergence rate. Fourth-order simulations yield nearly same accuracy (apart from minor gains in fault-slip and state misfits for well-resolved grids, h ≤ 30 m), regardless of whether t follows a constant CFL number (cfl = 0.5), or follows (33), with a CFL proportional to h [cfl = c * (h) < 0.5, for h < 100 m]. We conclude that MOSN's simple second-order treatment of time integration, with t set by the stability limit, does not require refinement. Fig. 6 further shows that a small amount of artificial viscosity (η = 0.025) marginally improves slip rate precision (convergence rate ∼1.4, versus ∼1.3 in the undamped case), but at the cost of significant degradation in the other variables, especially shear-traction and fault-slip (for which errors in the best resolved grids increase by ∼85 per cent and ∼300 per cent, respectively). We conclude that any gain from the suppression of spurious oscillations (most important in the slip-rate metric) is more than offset by reduced fidelity of peak tractions and accumulated slip, and in general we would not recommend use of artificial damping with MOSN, even in the second-order case.
Fourth-and mixed-order MOSN have accuracies very similar to that of BIEM for traction (apart from a single point at h = 20 m, which may be biased because the reference solution was from the BIEM family), state, and slip. However, in the better-resolved grids, BIEM performs systematically better than MOSN for slip rate (for example, fourth-order MOSN error is ∼57 per cent higher than BIEM error when h = 40 m). BIEM convergence rates are higher than both high-order MOSN implementations for slip rate and traction, but are still subquadratic. Convergence rates for slip and state are quadratic for both BIEM and MOSN.
We next study MOSN-Rosenbrock convergence rates on RS-A and RS-S problems. Fig. 7 summarizes results from the three spatial differencing schemes and excludes additional numerical devices (eqs 32 and 33). Since¯ differs between these test problems, we make the comparison only in terms of dimensionless resolution lengthN c, and omit h. In this case, errors for both slip rate (Fig. 7a) and traction (Fig. 7b) show very little dependence on state-evolution law (RS-A and RS-S errors differ by a few tens of percent), whereas both slip (Fig. 7c) and state-variable (Fig. 7d) errors split into separate populations for the respective evolution laws (RS-A and RS-S errors differ by roughly a factor 10). Convergence rates are given in Table 5, and show behaviour similar to that shown in Table 4 for problem RS-FH, namely In the RS-A test, BIEM misfits are represented by grey bullets and black-coloured circles, crosses, and triangles correspond to MOSN misfits from fourth-, mixed-and second-order simulations, respectively. Similar symbolism is followed to depict misfits in the RS-S case, but BIEM misfits are given by black symbols, while MOSN misfits are grey-coloured symbols. roughly quadratic convergence of all methods for slip and state variable, and mostly subquadratic but superlinear convergence for slip rate and traction.
We end this section by comparing accuracy and convergence rates of staggered MOSN-Euler and MOSN-trapezoidal methods with the corresponding MOSN-Rosenbrock results described above, limiting the comparison to the case of spatially fourth-order calculations and to problem RS-RH. Fig. 8 illustrates MAX-RMS misfits (of all fault variables), together with BIEM errors for reference, and Table 6 lists fitted convergence rates. We again use the numerical device of an artificially reduced Courant cfl number (eq. 33) with the MOSN-Euler scheme .  Figs 8(a) and (b) shows that all MOSN schemes display similar accuracy in modelling slip rate and shear traction, respectively, and Tables 4 and 6 indicate that they all display linear convergence with grid refinement. Figs 8(c) and (d) also show at most minor accuracy differences in MOSN solutions for fault slip and state variable, respectively, when h ≥ 40 m (N c ≤ 6). However, on finer grids time discretization errors become significant for MOSN-Euler, as shown by the fact that MOSN-Euler errors decrease when time step t follows a Courant number cfl < 0.5 [cfl = c * (h), given by eq. 33]. This device leads to a superquadratic convergence of the MOSN-Euler method that surpasses the quadratic decay of MOSN-trapezoidal and MOSN-Rosenbrock misfits and the superlinear convergence of MOSN-Euler with fixed cfl = 0.5 (see Tables 4 and 6). Thus, MOSN-Euler is the single case we have found in which time discretization errors can predominate, and then only under very limited conditions (slip and state variable only, and only for grids with very good spatial resolution).
Based on Figs 8(a)-(d), we conclude that the MOSN-trapezoidal scheme, with second-order convergence (for all fault variables) and accuracy comparable to MOSN-Rosenbrock, is our best candidate for practical implementations of R&S friction. MOSN-trapezoidal shows essentially the same error structure as MOSN-Rosenbrock, has substantial advantages over the latter in its simplicity of implementation, does not require substepping or interpolations, and also avoids the occasional artefacts of MOSN-Rosenbrock noted previously.

WAV E F O R M A S S E S S M E N T
We complete the study of MOSN precision by examining fault-plane waveforms. Zheng & Rice (1998) show that, for a background shear stress τ 0 below a critical value, τ pulse , mode III ruptures propagate as self-healing pulses rather than as expanding cracks, with τ pulse given by the minimum value of τ 0 such that a solution V exists for Above, f ss (V) is the steady state friction coefficient. We anticipate that mode II ruptures will show a transition to pulse-like behaviour in a manner similar to the mode III case. Test cases RS-A and RS-S have τ 0 > τ pulse according to eqs (8) and (34) and Table 1 (τ 0 = 60 MPa and τ pulse ∼ 55.9 MPa), so we anticipate crack-like propagation in those cases. To enable our test suite to sample pulse-like behaviour, RS-FH was designed such that τ pulse ∼ 28.312 MPa, exceeding the initial background stress τ 0 = 28 MPa (eqs 10, 34 and Table 1). Fig. 9 shows fault-plane time histories of slip, slip rate, shear traction and state variable for the three tests, at x = 6 km. The propagation is crack-like for both RS-A and RS-S: slip persists throughout the period of rupture expansion, and slip rate asymptotically approaches a steady-state level V ∞ (∼0.54 m s −1 for RS-A and ∼0.67 m s −1 for RS-S). Rupture propagation in RS-FH is pulse-like (slip ceases at a point while rupture is still expanding), showing that the mode III theory provides a good guide to the mode II transition to pulse-like behaviour as well. The residual stress levels for RS-A and RS-S correspond to the steady-state level of friction f SS at V ∞ given by (8), σ n f SS (V ∞ ) ∼ 55 MPa, (the difference in V ∞ between the two cases has no significant effect on f SS , due to the logarithmic dependence of f SS on V ). In case RS-FH, shear traction converges to the background level τ 0 = 28 MPa once it recovers from the abrupt drop associated with the breakdown process (Fig 9c). At fault location x = 6 km, the highest slip velocity peak occurs for model RS-FH, and the lowest for RS-A (Fig 9b). Of the three models, RS-FH induces the most abrupt decay of the state variable (Fig 9d), and, consequently, of traction (Fig 9c) at the rupture arrival, followed by RS-S, while RS-A shows the smoothest state and traction drops. In addition, Fig. 9 shows the effect of the state evolution model on the rupture speed, with RS-FH ruptures traveling the fastest and RS-A ruptures the slowest, with all tests having been done with similar input parameters (exceptions are indicated in Table 1). This state dependency of the rupture speed was previously depicted by Figs 3 and 4(a). Fig. 9 also shows the effects of grid resolution (measured byN c) on fourth-order MOSN time histories for each test case, and compares these with BIEM solutions of equivalent resolution (using MOSN-Rosenbrock waveforms as representative of fourth-order MOSN solutions for this qualitative analysis). For RS-A, we show the case h = 100 m (N c ∼ 11), for RS-S we show h = 40 m (N c ∼ 8), and for RS-FH we show h = 80 m (N c ∼ 3). The reference solution BIEM010 (N c ≥ 25, see Table 3), is also shown in each test case. Note that slip, traction and state waveforms, from both methods, BIEM and MOSN, are indistinguishable from the reference curve, even in the case with lowestN c, which is RS-FH. In particular, shear stress peaks at the S wave front and rupture front (Fig 9c), as well as the stress relaxation between these peaks, match the reference solution in all cases. At low resolution, as in the RS-FH example, both fourth-order MOSN and BIEM calculations of slip rate show spurious oscillations caused by grid dispersion as shown in Fig 9(b). Here, the inset illustrates that these oscillations are more pronounced for MOSN than for BIEM at the sameN c level. This behaviour is also characteristic of simulations with slip-weakening whenN c is of comparable magnitude (e.g. Day et al. 2005).

C O N C L U S I O N S
We have implemented rate-and state-dependent friction, including enhanced velocity-dependent weakening to emulate flash heating, in a split-node finite difference code (MOSN) with adjustable order of accuracy (fourth-, mixed-and second-order accuracy options). At a given split node, this formulation yields a coupled ODE system for slip velocity and state variable that becomes extremely stiff at low slip velocities (with stiffness inversely proportional to velocity squared, as shown in Appendix B). Numerical integration of the velocity-state system requires an implicit ODE solver, and the Rosenbrock method, having relatively high-order accuracy, provides useful benchmark solutions. That method is, however, sensitive to details of the elastic calculation scheme, and its adaptive substepping occasionally performs poorly in the presence of strong velocity weakening. In contrast, we obtain a robust, reliable and algorithmically simple solution scheme, free of substepping, and without loss of accuracy relative to the Rosenbrock scheme, by time-staggering slip velocity and the state variable, applying a low-order (trapezoidal rule) implicit ODE solver to the former, and performing semi-analytical (exponential) integration of the latter.
For a given error metric, finite difference solutions to all three rupture problems considered here (i.e. with aging, slip and flash heating state evolution laws) have similar convergence rates with respect to grid interval h. However, convergence rates differ among metrics associated with the fault-plane variables slip, slip rate, shear traction and state variable, respectively. Convergence rates for slip and state variable are roughly quadratic, while rates for slip velocity and traction are subquadratic. Convergence rates are similar for spatially fourth-, mixed-and second-order calculations, although the latter show some small differences from the other two, with a notable shift towards lower precision, especially when grid resolution is low (i.e. when the high stress-gradient zone at the rupture front, the cohesive zone, is resolved by fewer than about five grid nodes). Fourth-and mixed-order solutions have convergence rates for the state variable and the slip metrics that are very similar to the corresponding rates for BIEM solutions, and convergence rates somewhat lower than BIEM (by a factor of about 2/3) for the slip rate and traction metrics. The failure of (spatially) higher-order MOSN methods to significantly improve convergence rates relative to second-order methods is consistent with earlier results for slip-weakening friction models (Day et al. 2005;Rojas et al. 2008), so the limitation to second-order convergence rates found in the earlier studies is apparently not attributable simply to the presence of discontinuous derivatives in the shear-traction curves (since those curves have continuous derivatives in the friction model used in the current study). Nor are convergence rates limited by time-discretization errors, as long as the implicit velocity-state integration is at least second-order accurate (as the trapezoidal scheme is, but the backward Euler scheme is not), consistent with the second-order accuracy of the explicit integration used at interior nodes.
Self-healing pulse-like ruptures emerge in our mode II simulations under the same conditions as predicted by the mode III theoretical results of Zheng & Rice (1998), i.e. when initial shear stress is lower than a critical value (τ pulse ) given by the zero-velocity intercept of the radiation damping line tangent to the steady-state velocity-weakening curve. Fault-plane time histories (slip, slip velocity, shear traction and state variable) from MOSN simulations of both pulse-like and crack-like ruptures are nearly indistinguishable from the corresponding high-resolution BIEM reference solutions, apart from small, spurious oscillations in slip velocity for the lowest-resolution cases (with just three-node cohesive-zone resolution). These slip-velocity oscillations are less pronounced for rate-and-state problems than they are for slipweakening problems (given similar grid resolution), because of the natural damping provided by the direct velocity-dependent term in the former. As a result, artificial viscosity may not be advantageous in simulations with rate and state friction, even for the lowest (second) order schemes, contrary to our conclusions (e.g. Day et al. 2005; for low-order simulations with slip weakening friction.

A C K N O W L E D G M E N T S
This work was supported by the National Science Foundation, under grant EAR-0810271 and EAR-0623704, and by the Southern California Earthquake Center (SCEC). SCEC is funded by NSF Cooperative Agreement EAR-0529922 and USGS Cooperative Agreement 07HQAG0008. Otilio Rojas was also funded by the Computational Science Research Center (SDSU) and the Universidad Central de Venezuela (UCV-CDCH).
(ii) Fourth-order case: (iii) Mixed-order case: In the case (ii), we have corrected typographical errors in the components of fourth-order operators G (K) and D (K) that appeared in Rojas et al. (2008). Note that in the case (iii), the order of accuracy of approximations Gf and Dv reduces in the vicinity of nodes z −M , z 0 and z N .

A P P E N D I X B : S T I F F N E S S A N D T H E V E L O C I T Y-S TAT E S Y S T E M
In this appendix, we discuss and justify the time-integration procedures, via a linear stability analysis of the spatially discretized governing equations, that we use to couple rate-and-state friction laws with an elastic medium. The essential ideas are captured in a 1-D setting in which fields in the medium vary only in the fault-normal direction, slip is unidirectional, and symmetry conditions are assumed about the fault so we can restrict attention to one side only.
We consider general expressions for the friction coefficient, f (V , ψ), and evolution function, G(V , ψ). The exterior boundary condition is one of constant stress, τ 0 , though this is purely for later convenience and the general results of this analysis hold for more general boundary conditions at z = W . The initial conditions are consistent with the boundary conditions, that is, τ 0 = σ n f (V 0 , ψ 0 ). To simplify later expressions, let v x (z, t) = v 0 + v (z, t), τ xz (z, t) = τ 0 + σ (z, t) and ψ(t) = ψ 0 + (t). Note that v(z, t), σ (z, t) and (t) satisfy with homogeneous initial conditions v (z, 0) = 0, σ (z, 0) = 0, and (0) = 0, and boundary conditions We spatially discretize by dividing the domain 0 ≤ z ≤ W into N cells of width h = W /N . Velocity is defined at the N + 1 cell edges: Additionally, we define two extra values of stress at z = 0 and z = W , which are denoted σ 0 (t) and σ N+1 (t), respectively, which are used to enforce the boundary conditions by setting Spatial derivatives of stress at the velocity gridpoints are approximated by the (N + 1) Likewise, spatial derivatives of velocity at the interior stress gridpoints are approximated by the N × The spatially discretized governing equations are obtained by first replacing σ 0 (t) and σ N+1 (t) in (B11) using (B9) and (B10), and then substituting (B11) and (B12) into (B4). The resulting set of equations, together with the discretized version of (B5), are These equations constitute a non-linear system of ODEs. To proceed further in the analysis, we must linearize the two non-linear functions f (V , ψ) and G(V , ψ) by assuming |v 0 | |V 0 | and | | |ψ 0 |. Expanding to first order gives Rupture modelling with strong velocity-weakening friction 1855 where the partial derivative terms are evaluated at V = V 0 and ψ = ψ 0 . Substituting these expressions into (B13)-(B15) yields a linear system of ODEs The stability of the system of ODEs (B18) is assessed by examining the 2N + 2 eigenvalues, λ i , of the matrix [M ij ]. An instructive starting point for discussing the eigenvalues of [M ij ] is the case of a fault that slides at a constant coefficient of friction with no state evolution, that is, ∂ f /∂ V = ∂ f /∂ψ = ∂G/∂ V = ∂G/∂ψ = 0. For this special case, and for the mimetic differentiation matrices [G ij ] and [D ij ], all eigenvalues lie on the imaginary λ-axis. The magnitudes of the largest eigenvalues determine the maximum permissible time step when using an explicit time-stepping method. For the leapfrog method as applied in this work, the maximum CFL ratio, expressed in terms of the shear-wave transit time across a cell h/c s , is 2c s /(h max i |λ i |). In this 1-D geometry, the maximum CFL ratios for the second-order, mixed-order and fourth-order differentiation matrices are 0.9306, 0.8526 and 0.8165, respectively.
When we have non-trivial friction laws and evolution equations, the eigenvalues are shifted relative to their values for the homogeneous boundary conditions. To aid the following analysis, the partial derivatives appearing in (B18) are ∂ f /∂ V = a/V, ∂f /∂ψ = 1, in which and We illustrate the main features of (B18) by focusing on the specific case of RS-S. The eigenvalues of [M ij ] are shown in Fig. B1 for several values of V and f (V , ψ). For rate-and-state friction laws, most the eigenvalues are shifted slightly off of the imaginary λ-axis into the left half-plane. While this corresponds to a slight damping of the eigenvectors, it does imply a very mild numerical instability when time stepping is performed with the leapfrog method. The leapfrog (or implicit midpoint) method is only stable when used to solve ODE systems having imaginary eigenvalues since its absolute stability region consists only of a finite segment of the imaginary axis centered on the origin [e.g.
An ODE system is said to become stiff when different components evolve over different timescales, and an appropriate way to quantify this is to define stiffness as S = |λ − |/ |λ + |. Stiffness occurs in our problem when t a t ψ or, equivalently, 1/ξ 1, which occurs when slip velocity is sufficiently small, i.e. V √ κaσ n L/δ ≈ 1m/s. We can expand (B27) and (B28) in powers of the small parameter 1/ξ to obtain, to leading order, the following approximations: The eigenvalue associated with rapid damping of velocity perturbations by the direct effect, λ − , is identical, in this limit, to the very large negative eigenvalue in the spatially discretized 1-D problem discussed at the beginning of this appendix and shown in Fig. B1(a). The other eigenvalue, λ + , can be either negative or positive; the sign is determined by that of η/δ−1. For perturbations about steady sliding (at speeds below V w for RS-FH), η/δ − 1 = b/a − 1; the sign of λ + is then determined by whether the steady state friction coefficient is velocity-weakening or velocity-strengthening. We conclude from our analysis that except during conditions of rapid coseismic sliding, the two eigenvalues differ by many orders of magnitude and the system is very stiff. Stiffness is highly sensitive to slip velocity (S ∝ V −2 ), and exhibits considerable variation during the course of a single simulation (ranging from S ∼ 10 25 at V ∼ 10 −12 m s −1 to S ∼ 1 at V ∼ 1 m s −1 ).
The difficulty with stiff systems is that explicit numerical integration is only stable if transient evolution over all involved timescales is resolved (i.e. for the velocity-state system, time steps must be smaller than both t a and t ψ ). The equations are stiff when t a t ψ , but we are only interested in the system over the longer time scale t ψ , which is comparable to or less than the time step required for stable and accurate integration of wave propagation in the interior. Thus we must use an implicit numerical method to solve the velocity-state system to avoid the time step restrictions imposed by explicit methods. We have successfully implemented a variety of these methods, ranging from high-order Rosenbrock methods to simpler methods like implicit Euler and the implicit trapezoidal rule; these are discussed in the text.
In Fig. B2, we show how the two eigenvalues evolve as the rupture passes a point on the fault. This is shown using the time histories of fields in the BIEM010 solutions to the RS-A and RS-S test problems at fault location x = 12 km; results are similar at other points  Figure B2. Evolution of the eigenvalues of the velocity-state system, as given by eq. (B27), as a function of slip rate V at x = 12 km for the BIEM010 solutions to test problems RS-A and RS-S. (a) One eigenvalue, λ + , is associated with evolution over a time scale longer than the time step required for stable and accurate integration of wave propagation ∼h/c s . (b) The other eigenvalue, λ − , is associated, in the small V limit, with the rapid damping of velocity perturbations by direct effect. Stable integration with time steps ∼h/c s necessitates the use of methods designed for stiff systems. Note log-log scale. −6.72 × 10 11 , −3.09 × 10 −13 ) ( −6.72 × 10 7 , −2.02 × 10 −9 ) ( −6.72 × 10 3 , −9.45 × 10 −6 ) ( −6.72 × 10 −1 , −2.48 × 10 −3 ) ( −7.96 × 10 −2 ± 3.37 × 10 −1 i −3.40 × 10 −1 ± 2.62 × 10 −1 i −6.72 × 10 2 , −2.37 × 10 −3 ) (  Figure B3. Evolution of velocity and friction coefficient at x = 12 km for the BIEM010 solution to test problem RS-FH, with eigenvalues given at several times (labelled A-G). As V increases from the initially low background value and frictional resistance increases by the direct effect (A-D), the two eigenvalues approach each other in magnitude. As the rupture front passes and state evolution occurs, the friction coefficient drops (E-F). During this time, the two eigenvalues are complex conjugates with negative real part. After the slip pulse passes (G), the eigenvalues return to real values.