Abstract. We present an efficient and rigorous numerical procedure for calculating the elastodynamic response of a fault subjected to slow tectonic loading processes of long duration within which them are episodes of rapid earthquake failure. This is done for a general class of rate-and state-dependent friction laws with positive direct velocity effect. The algorithm allows us to treat accurately, within a single computational procedure, loading intervals of thousands of years and to calculate, for each earthquake episode, initially aseismic accelerating slip prior to dynamic rupture, the rupture propagation itself, rapid post seismic deformation which follows, and also ongoing creep slippage throughout the loading period in velocity-strengthening fault regions. The methodology is presented using the two-dimensional (2-D) antiplane spectral formulation and can be readily extended to the 2-D in-plane and 3-D spectral formulations and, with certain modifications, to the space-time boundary integral formulations as well as to their discretized development using finite difference or finite element methods. The methodology can be used to address a number of important issues, such as fault operation under low overall stress, interaction of dynamic rupture propagation with pore pressure development, patterns of rupture propagation in events nucleated naturally as a part of a sequence, the earthquake nucleation process, earthquake sequences on faults with heterogeneous frictional properties and/or normal stress, and others. The procedure is illustrated for a 2-D crustal strike-slip fault model with depth-variable properties. For lower values of the state-evolution distance of the friction law, small events appear. The nucleation phases of the small and large events are very similar, suggesting that the size of an event is determined by the conditions on the fault segments the event is propagating into rather than by the nucleation process itself. We demonstrate the importance of incorporating slow tectonic loading with elastodynamics by evaluating two simplified approaches, one with the slow tectonic loading but no wave effects and the other with all dynamic effects included but much higher loading rate.

We study the stability of steady sliding between elastically deformable continua using rate and state dependent friction laws. That is done for both elastically identical and elastically dissimilar solids. The focus is on linearized response to perturbations of steady-state sliding, and on studying how the positive direct e ect (instantaneous increase or decrease of shear strength in response to a respective instantaneous increase or decrease of slip rate) of those laws allows the existence of a quasi-static range of response to perturbations at su ciently low slip rate. We discuss the physical basis of rate and state laws, including the likely basis for the direct e ect in thermally activated processes allowing creep slippage at asperity contacts, and estimate activation parameters for quartzite and granite. Also, a class of rate and state laws suitable for variable normal stress is presented. As part of the work, we show that compromises from the rate and state framework for describing velocity-weakening friction lead to paradoxical results, like supersonic propagation of slip perturbations, or to ill-posedness, when applied to sliding between elastically deformable solids. The case of sliding between elastically dissimilar solids has the inherently destabilizing feature that spatially inhomogeneous slip leads to an alteration of normal stress, hence of frictional resistance. We show that the rate and state friction laws nevertheless lead to stability of response to su ciently short wavelength perturbations, at very slow slip rates. Further, for slow sliding between dissimilar solids, we show that there is a critical amplitude of velocity-strengthening above which there is stability to perturbations of all wavelengths. ?

[1] The spontaneously propagating shear crack on a frictional interface has proven to be a useful idealization of a natural earthquake. The corresponding boundary value problems are nonlinear and usually require computationally intensive numerical methods for their solution. Assessing the convergence and accuracy of the numerical methods is challenging, as we lack appropriate analytical solutions for comparison. As a complement to other methods of assessment, we compare solutions obtained by two independent numerical methods, a finite difference method and a boundary integral (BI) method. The finite difference implementation, called DFM, uses a traction-at-split-node formulation of the fault discontinuity. The BI implementation employs spectral representation of the stress transfer functional. The three-dimensional (3-D) test problem involves spontaneous rupture spreading on a planar interface governed by linear slip-weakening friction that essentially defines a cohesive law. To get a priori understanding of the spatial resolution that would be required in this and similar problems, we review and combine some simple estimates of the cohesive zone sizes which correspond quite well to the sizes observed in simulations. We have assessed agreement between the methods in terms of the RMS differences in rupture time, final slip, and peak slip rate and related these to median and minimum measures of the cohesive zone resolution observed in the numerical solutions. The BI and DFM methods give virtually indistinguishable solutions to the 3-D spontaneous rupture test problem when their grid spacing Dx is small enough so that the solutions adequately resolve the cohesive zone, with at least three points for BI and at least five node points for DFM. Furthermore, grid-dependent differences in the results, for each of the two methods taken separately, decay as a power law in Dx, with the same convergence rate for each method, the calculations apparently converging to a common, grid interval invariant solution. This result provides strong evidence for the accuracy of both methods. In addition, the specific solution presented here, by virtue of being demonstrably grid-independent and consistent between two very different numerical methods, may prove useful for testing new numerical methods for spontaneous rupture problems.Citation: Day, S. M., L. A. Dalguer, N. Lapusta, and Y. Liu (2005), Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture,

No abstract

1 Supplementary Information S1. Model parameters Figure S1 and Table S1 show the model geometry and parameters used in the simulation presented in the main article. S2. Model response that reproduces a range of observationsThe rich model response described in the main text is animated by the two supplementary movies which show the time evolution of slip rate distribution V z (x, z, t) on the fault from 4300 to 4800 years (movie I) and from 3500 to 4000 years (movie II) after the beginning of the numerical simulation. Coseismic slip rates of the order of 0.1 m/s and higher are indicated by white to blue colors. The slower slip rates are shown in progressively darker yellow and orange, with the long-term plate rate of 10 -9 m/s indicated in red. Locked regions which have, by definition, much smaller slip rates than the plate rate are indicated by black patches.The simulated time and the maximum slip rate over the fault at that time are given at upper-right and upper-left corners of the movies, respectively. The movie frames are taken every 30 computational time steps, which vary in the computation, and hence the time intervals between the frames are quite heterogeneous. During coseismic periods, the typical time step is 2/150 s and the typical interframe rate is 0.4 s. In the interseismic periods, the interframe rates reach values of about 1 year.Movie I illustrates the behavior of the model that combines, in the same patch, creep in the interseismic periods and large seismic slip. It starts at a time when the interseismic fault slip rate distribution is as expected for the slow-rate friction properties adopted in the model: patch A is mostly locked and patch B is creeping with the long-term plate rate, like the creeping areas that surround the two patches. During the following interseismic period, the locked region of patch A shrinks by penetration of creep from the surrounding areas, which occurs due to stress concentration at the boundary between creeping and locked regions (e.g., Tse and Rice, 1986;Lapusta et al., 2000). Earthquake rupture nucleates when the creeping region within patch A becomes comparable to the nucleation size estimates (Rice and Ruina, 1983;Rice et al., 2001;Rubin and