and Tang (11) on a simple cellular-automata model of failure at a critical stress and stress shedding to neighboring cells.

Rice (5) suggested that some of the above models [not including the work of Horowitz and Ruina (9)] showed aspects of slip complexity as an artifact of their use of simplified constitutive laws, which are not fully compatible with continuum elastic models, as we discuss below. It was also suggested that such models produce results that depend on the cell size (or in the B-K context block spacing) *h* of the numerical scheme and should not properly be regarded as models of smooth faults. Instead, they may be taken as approximate models of strongly disordered faults with a fundamental scale of heterogeneity equal to *h*. Such an interpretation has been expanded upon by Ben-Zion and Rice (6, 12) whose results are summarized subsequently. The conclusion reached was that some of the modeling introduced as showing complexity through the first route, dynamics on a smooth fault, may have a surer physical interpretation as examples of the second route to complexity, strong heterogeneity. This conclusion has, however, been contested (13, 14) and we return to it below.

It has been a concern in evaluating the first possibility, *(i)* that previous studies (5, 6) from this group used a quasidynamic approximation, that we explain shortly, rather than fully solving the equations of elastodynamics, and that earlier studies (4) used similar quasistatic analyses. The question was whether the conclusions (5, 6, 12) reached on complexity would change if there was full inclusion of inertial dynamics. Indeed, Carlson and Langer (10) suggested that the complexity of slip history predicted for their B-K array was linked to inertial dynamic effects, which were fully included (within the lumped mass and one-dimensional fault model approximations) in their simulations. As we remark below, however, we think there is evidence that features of their classical constitutive description, and also of their numerical solution algorithms [abrupt strength drop “*σ*” imposed at onset of slip, a feature repeated in Shaw (13) and Myers *et al.* (14)], which are not fully suitable for representation of slip on a smooth fault, may have a dominant role in giving the small-event complexity of the G-R type that they found. Cochard and Madariaga (15, 16) have also argued, from elastodynamic studies of repeated antiplane slip ruptures on an isolated fault segment in an elastic continuum, that inertial effects are critical to the complexity they observe. However, they seem to require a restricted class of constitutive assumptions to find such complexity and the small event distribution they show is not of G-R type (16).

We have recently developed the methodology to use full elastodynamics in some of our modeling of repeated earthquakes and present some results that support previous conclusions (5, 6, 12) from this laboratory; we have found no evidence that inclusion of fully inertial elastodynamics brings small event complexity of G-R type to our smooth fault models.

To understand the response of a smooth fault, we have studied slip-rupture sequences along a smooth vertical strike-slip fault in an elastic half-space. The basic model for this, from Rice (5) and subsequent studies (6, 17, 18), is shown in Fig. 1. The fault is driven below a depth of 24 km at rate *V*_{p1}(=35 mm/year in cases shown here), and the slip history is calculated in the shallower zone. The modeling includes laboratory-constrained temperature, and hence depth, variability of friction properties, based on experiments on granite gouge under hydrothermal conditions by Blanpied *et al.* (19). The laboratory data were adapted (5) for a depth-variable model of frictional properties of the crust by using a San Andreas fault geotherm, thus making the friction properties, denoted *a* and *b* subsequently, functions of depth *z* in the crust. The depth variation of the friction properties implies a transition from potentially

unstable velocity-weakening behavior, *b>a,* to inherently stable velocity strengthening, *b<a,* at an ≈14-km depth, and this serves to confine unstable slip to a roughly 15-km seismogenic depth. Over the 4.0- to 13.5-km depth, *a*=0.015 and *b*=0.019 in the model (5).

The friction model incorporates the observed displacement-dependent features of friction as in rate- and state-dependent frictional constitutive laws. Such assures that numerical fault models have a well-defined limit with refinement of the computational grid (5); the same is not necessarily the case with modeling based on classical friction laws without gradual displacement-dependent features of strength transition. In the two-dimensional (2D) versions of the modeling, slip *δ* varies just with depth *z, δ*=*δ(z, t).* In three-dimensional (3D) versions, it varies with depth and distance *x* along strike, *δ*= *δ(x, z, t),* with the slip distribution being repeated periodically along strike with wavelengths ranging from 240 to 960 km. Here *δ(x, z, t)=u*_{x}(*x,* 0^{+}, *z, t*)*−u*_{x}(*x,* 0^{−}, *z, t*)*,* where *(u*_{x}*, u*_{y}*, u*_{z}*)* is the displacement vector and we disallow discontinuities in *u*_{y} and *u*_{z}*.*

Our procedure for fully elastodynamic analysis is based on the spectral method of Perrin *et al.* (20), who use it for 2D studies of spontaneous antiplane rupture with rate- and state-dependent constitutive laws. The formulation has been extended (21) to spontaneous rupture in 3D, with slip or crack opening varying with position *x, z* on a plane and with time, but implemented thus far only for problems of tensile cracks. Thus all of our fully dynamic earthquake simulations to date are of 2D type, *δ*=*δ(z, t).*

The shear stress *τ(x, z, t),* corresponding to the tensor stress component *τ*_{yx}(*x,* 0, *z, t*) on the fault, has an elastodynamic representation in the form

*τ(x, z, t)=τ°(x, z)*+*ϕ*[*δ(x′, z′, t′)−V*_{p1}*t′; x, z, t*] −*μV(x, z, t)*/2*c*_{s},

where *τ*^{o}*(x, z)* is an initial value, *μ* is shear modulus, *V=∂δ/∂t, c*_{s} is shear wave speed, and *ϕ* represents a quantity at *x, z, t,* which is a linear functional of its first argument, *δ(x', z', t')−*V_{p1}*t′*, over all *x′, z′, t′* values within the “wave cone” influencing *x, z, t*. The functional *ϕ* can be written as an integral over that wave cone, as in several sources on elastodynamics (22, 23, 24, 25), and such is the usual starting point for related boundary integral formulations. We represent both *ϕ* and *δ*−*V*_{p1}*t,* in the form of Fourier series in *x* and *z* and use expressions from elastodynamics for how the time-dependent Fourier coefficients for *ϕ* are given as convolutions over prior time of those for *δ*−*V*_{p1}*t* .

The route to handling entire earthquake cycles elastodynamically begins with the recognition (20) that