stably creeping lower region of the fault. The period of simple harmonic motion associated with this term is proportional to the time required for an elastic wave to traverse the crust depth, which is the same as the characteristic duration of slip at a single point on the fault in a very large, characteristic event. For real faults, this period is the order of seconds, the wave speed is the order of 10 kilometers/sec, and the corresponding length scale is the order of 10 kilometers. Accordingly, our unit of position x is this seismogenic depth.
The last term in Eq. 1 is the velocity-weakening stick-slip friction, which we usually have taken to have the form
(For some applications, we have found piecewise linear or exponential forms of the velocity-weakening function to be more convenient. The qualitative properties of the model seem to be insensitive to the specific functional form of Φ.) We have chosen the slipping threshold to be unity by scaling U so that this threshold is reached in a completely uniform system when the displacement U is equal to −1 for all x. Thus our displacements U are measured in units of the characteristic slip in a very large event, a length of order meters for real faults. With this scaling, the dimensionless loading rate v is the ratio of the true loading rate, centimeters per year, to the characteristic slipping speed, meters per second. Therefore, v≈10−8.
Note that Eq. 2 is a specially simple friction law with a sharply defined slipping threshold and no extra rate or state dependencies (7, 8). This friction law permits no stable creep at small velocities, and the slipping threshold is fully restored immediately upon resticking.
The only important system-dependent parameter in this model is the velocity-weakening rate αv, which cannot be measured directly in the Earth. In the absence of information to the contrary, we assume that αv must have a value of order unity. The dynamic instability associated with the velocity-weakening force law is the principal cause of slip complexity in this system. A larger value of αv means less energy dissipation and stronger instability; smaller αv means more dissipation and weaker instability. We can see this explicitly from a simple linear analysis of the first-order displacement δU(x, t) during a slipping motion. Let δU~exp(ikx+ωt), where k is a wavenumber and ω is the amplification rate. Then, from Eqs. 1 and 2 we find that
This instability persists—ω remains positive—down to arbitrarily small wavelengths. As a result, the model needs some sort of short-wavelength cutoff. We shall return to this point.
The parameter σ that appears in Eq. 2 is the force drop or, equivalently, the acceleration of a block at the instant when slipping begins. We have introduced α in this version of the model primarily as a technical device that allows us to study the limit v → 0. In the absence of σ, events begin infinitely slowly in that limit, which obviously is inconvenient for numerical purposes. With a value of σ of order 10−2 and vanishingly small v, on the other hand, we need not carry out explicit time integrations between slipping events and thus can study the system for very large numbers of loading cycles. We shall see, however, that σ has greater physical significance in the two-dimensional slip-weakening model.
We have performed numerical integrations of the equation of motion Eq. 1 over the equivalent of many thousands of seismic loading cycles with systems containing thousands of blocks. Our most important observation is that, in this limit in which the loading rate v is very small, the system exhibits persistent slip complexity. That is, if we start with any slightly nonuniform configuration U(x, t=0) and integrate Eq. 1 forward in time, without introducing any further noise or irregularities, we find chaotic motion. A portion of one such simulation is illustrated in Fig. 1, where we plot a sequence of fully stuck configurations. The displacements U(x) jump forward during slipping events, and the areas swept out by these events are our dimensionless seismic moments, M=∫δUdx.
Our extensive numerical studies indicate that what we are seeing in Fig. 1 is a finite sequence of configurations that belongs to a statistically stationary distribution of event sizes. Stationarity means that, after initial transients have disappeared, this distribution is entirely independent of the initial configuration. For αv≥2, the distribution contains two distinct populations: a region of small to moderately large localized events (corresponding to the smaller events that cluster in the local minima in Fig. 1) that obey something like a Gutenberg-Richter (GR) law; and a region of large delocalized events—irregularly recurring characteristic earthquakes—whose frequency exceeds that of the extapolated GR law, and that account for essentially all of the moment release.
A typical frequency-magnitude distribution is shown in Fig. 2. Here, R(μ) is the differential distribution of frequencies of events of magnitude μ=lnM. The distribution of localized events in this log-log plot has slope −1, a value which persists for all αv>2. For smaller αv, we find smaller slopes; and we also find that, although complexity persists, the characteristic-event peak disappears near αv=1 as the instability becomes weaker.
The large, delocalized events consist of pairs of slipping pulses that emerge from a nucleation zone and propagate in opposite directions along the fault at roughly the elastic wave speed. We have studied these Heaton pulses (9) in considerable analytic detail (10, 11) and find that they behave in many respects like shock fronts whose widths are controlled by the short-wavelength cutoff. The anomalous sharpness of these pulses may be an important factor in generating the small-event complexity that is seen in these models. The crossover between localized and delocalized events occurs at an event size that is proportional to the crust depth multiplied by an amplification factor (order of 2 or 3 for the parameters used