# Positive feedback, memory, and the predictability of earthquakes

^{*}Department of Earth Sciences, University of Southern California, Los Angeles, CA 90089-0740;^{‡}Institute of Geophysics and Planetary Physics, and Department of Earth and Space Science, University of California, Los Angeles, CA 90095; and^{§}Laboratoire de Physique de la Matière Condensée, Centre National de la Recherche Scientifique Unité Mixte de Recherche 6622, and Université des Sciences, B.P. 70, Parc Valrose, 06108 Nice Cedex 2, France

See allHide authors and affiliations

## Abstract

We review the “critical point” concept for large earthquakes and enlarge it in the framework of so-called “finite-time singularities.” The singular behavior associated with accelerated seismic release is shown to result from a positive feedback of the seismic activity on its release rate. The most important mechanisms for such positive feedback are presented. We solve analytically a simple model of geometrical positive feedback in which the stress shadow cast by the last large earthquake is progressively fragmented by the increasing tectonic stress.

Are earthquakes predictable? The answer, of course, depends on what is meant by a prediction. In the broadest sense, the plate tectonics paradigm makes predictions. It predicts that earthquakes are far more likely to occur at the boundaries between plates than within their interiors. Actually, plate tectonics theory was in part based on this “in-sample” observation, which is verified continuously “out-of-sample.” It also predicts an overall rate to the process. Averaged over time, the summed moments of the earthquakes is consistent with the relative motion between the plates determined from the analysis of magnetic anomalies, correcting for aseismic visco-plastic deformations.

The forecasting of individual large events has been more problematical. Although the paleoseismological dating of large prehistoric earthquakes has confirmed the plate tectonics hypothesis, the timing between individual events is extremely erratic. For example, trenching by Sieh *et al.* (1) found that the average recurrence interval for the last ten large earthquakes on the San Andreas Fault north of Los Angeles is about 132 years. Because the long-term slip rate on the Southern San Andreas is about 3 cm/year, this corresponds to an average displacement of about 4 m per large earthquake—a very reasonable value. The problem is that the intervals between events range from 44 to 332 years. This lack of quasiperiodicity in large events is also evident in the failed Parkfield prediction, where the current interval is approaching 35 years despite a prior string of intervals averaging about 21 years (see, e.g., ref. 2). These observations have dimmed the hope that large earthquakes can be forecast based solely on the past history of large events on the same fault.

An alternative forecasting strategy is based on physical precursors observed to occur just before macroscopic failure in the laboratory. Most of these precursors are associated with microfracture damage and the associated dilatancy observed to precede the formation of a macroscopic shear failure of rock specimens under compressive loading. These laboratory observations have been incorporated into the “dilatancy-diffusion model” for earthquakes (3–6). However, the search for physical precursors before large earthquakes has been disappointing. The high hopes raised by the reports of Chinese success in using physical precursors to forecast the 1975 Haicheng earthquake have dissipated with the worldwide failure to produce additional valid predictions (for a review, see ref. 7).

In the United States, current work on earthquake prediction is primarily based on the search for precursors to large events in the seismicity itself. One motivation comes from a statistical physics interpretation of regional seismicity as being characteristic of a system at or near a statistically stationary dynamical critical point dubbed self-organized criticality (see ref. 8 for a review of mechanisms). Such a self-organized critical state is characterized by power law distributions of event sizes and long-range spatial correlation (9) of fluctuations around the statistically stationary state. Because earthquakes are characterized by several power laws, such as the Gutenberg–Richter frequency–magnitude relation, Omori's Law for aftershocks, and the fractal distribution of hypocenters, the application of this concept of self-organized criticality to earthquakes is now often taken for granted in the seismological community.

However, the implication of this self-organized critical (statistically stationary) state for the predictability of large earthquakes remains controversial. Geller *et al.* (10) pointed out that if the crust follows the strict tenets of self-organized criticality, then earthquakes are inherently unpredictable. According to this view point, simple self-organized critical systems having short-range interactions and conserved stress (except at the boundaries) tend to sit at or near the statistically stationary dynamical critical point such that the probability for a small event to cascade into a large one is not preconditioned by previous events. Physically speaking, this is because the long-range stress correlations are not significantly destroyed by large events in such systems. This view has actually been proved wrong by direct tests of the predictability of large events in several sandpile models, which show detectable correlations between precursory smaller events and the largest events (11), though the correlations may be very weak in some cases. The physical reason lies in the fact that a statistically stationary dynamical critical state is in fact the result of subtle and long-range correlations in space and time between events. As a consequence, a large event occurs on a long-range correlated landscape. In practice however, these effects are small and hard to detect and quantify precisely.

The alternative point of view is that large system-wide events produce significant perturbations that move a system away from the statistically stationary dynamical critical point. Nonconservative automata (see, e.g., ref. 12) or automata having correlated heterogeneity on the scale of the network have been shown to produce such perturbations (13, 14). Because these systems have a strong memory of past events, prediction becomes much more possible. Pushing this concept, Sornette and Sammis (15) have proposed that great earthquakes are themselves “critical points”—i.e., the culmination of a nonstationary precursory activity accelerating up to the critical point. We emphasize that this concept is fundamentally different from the long-time view of the crust as evolving spontaneously in a statistically stationary critical state, called self-organized criticality (SOC). In the SOC view, all events belong to the same global population and participate in shaping the self-organized critical state. By contrast, in the critical point view a great earthquake plays a very special role and signals the end of a cycle on its fault network. The dynamical organization is not statistically stationary but evolves as the great earthquake becomes more probable. As we said, predictability becomes possible and may be implemented by monitoring the approach of the fault network toward a critical state.

The basic tenets of the critical point model for regional seismicity can be simply stated as follows: (*i*) A large (network-sized) earthquake is only possible when the crust has reached a critical state. Highly stressed patches must be correlated at the scale of the fault network if an event, once nucleated, is to grow by jumping geometrical and rheological barriers. (*ii*) A large earthquake moves portions of the network out of the critical state by destroying stress correlations. This lack of correlation produces the observed shadow zones for intermediate size events. (*iii*) Tectonic loading combined with the stress transfer from smaller events reestablish long-range stress correlations, thus making the next large event possible. Note that nothing in this hypothesis implies that the seismic cycle, thus defined, has to be periodic. Presumably, observed variations in recurrence intervals discussed above result from details of the stress transfer from smaller events and from details of the nucleation process. Also, nothing in the hypothesis says that the large earthquake has to occur when the system reaches the critical state—only that a large earthquake is then possible. The exact time of the large earthquake may depend on details of the nucleation process, which may itself have a significant time dependence (see, e.g., ref. 16). Predictability in this model comes from monitoring the approach of the fault network back to the critical state.

Observations of the acceleration of seismic moment leading up to large events and “stress shadows” following them have been interpreted as evidence that seismic cycles represent the approach and retreat of a fault network for a critical state. Harris and Simpson (17) and Jones and Hauksson (18) present observational evidence for “Coulomb-stress shadows” following the 1857 Fort Tejon, the 1952 Kern County, and the 1992 Landers earthquakes. Harris and Simpson found *no* M > 5.5 earthquakes in the Coulomb-stress shadows for 50 years following the 1857 and 1952 events, whereas M > 5.5 events continued at a normal rate in areas where the Coulomb stress was not lowered. This lack of intermediate-sized events is the expected result of a reduction in the stress correlation length associated with a retreat from the critical state following a large event.

The fundamental hypothesis of the critical point model is that fault networks naturally evolve back toward the critical state following a large event, and that this approach to criticality can be monitored by using regional seismicity. Three approaches to this monitoring have been proposed. The first is to monitor the increase in event size as a proxy for the growing correlation length. Sykes and Jaumé (19) documented this effect as an increase in intermediate-sized events in a large region surrounding the impending large earthquake. Bufe and Varnes (20) quantified the effect of these intermediate-sized events in terms of a power law time-to-failure equation of the form where *E*_{j} is the energy of the *j*th event, and *t*_{c} is the time that the system reaches the critical state and a large event becomes possible. The quantity *∑**E* is usually referred to as the “cumulative Benioff strain.” The parameter *B* is negative and *m < 1* (typically *m ≈ 0.3)* so that the cumulative Benioff strain accelerates up to the finite value *A* with a diverging slope as *t* goes to *t*_{c}. In practice, only the intermediate-sized events (typically two magnitude units below the large event or greater) are included in the sum **1** because the Benioff strain is dominated by small events.

Sornette and Sammis (15) pointed out that Eq. 1 is the expected behavior near a critical point and proposed in addition to enrich it with log-periodic corrections (see also ref. 21). Bowman *et al.* (22) optimized the fit to Eq. 1 with respect to the size of the critical region to find a simple scaling relation between magnitude and radius for circular regions. In fact, they found the same scaling proposed by Keilis-Borok and Malinovskaya (23) based simply on the width of the fault network. A good review of this line of research is in ref. 24.

More recently, a second method has been proposed by Zöller *et al.* (25) who estimate the growing correlation length directly from the earthquake catalog by using single-link cluster analysis. Their method has the advantage of working with smaller events, thus improving the statistical significance of the results.

A third method for monitoring the approach to criticality is based on the notion that a system approaching the critical state becomes increasingly sensitive to tidal stress perturbations. Yin *et al.* (26) proposed that a correlation between solid earth tides and seismicity develops locally just before a large earthquake. They quantify this effect by calculating the average seismic energy release in a region surrounding the impending earthquake during periods when the effective tidal stress on the fault plane is positive divided by the energy released during periods when it is negative. Dubbed the Load/Unload Response Ratio (LURR), Yin *et al.* (26, 27) reported LURR calculations for a number of earthquakes in China, Australia, and California, and suggest that this ratio increases prior to the large earthquakes that they studied. The time period over which this increase occurs appears to be several months to 1 or 2 years. Recently, Yin *et al.* (27) have broadened their rationale to encompass the critical point model for large earthquakes. They support this interpretation with examples that illustrate that the largest and best-defined anomalies in LURR seem to occur when the choice of the region over which seismicity is averaged follows the scaling law developed by Bowman *et al.* (22), using the power law acceleration of regional seismicity. In a related development, Wang *et al.* (28) carried out experiments with the “lattice solid model” of Mora and Place (29, 30) in which sinusoidal perturbations to the loading force were applied in their numerical simulations. The results showed a confirmation that anomalous values of LURR occurred prior to large events in the model.

There are, however, theoretical arguments that the LURR effect should not be observable. Dieterich (16, 31) points out that fault patches with rate- and state-dependent friction properties are characterized by an interval of self-driven accelerating slip prior to instability. This delayed nucleation provides an explanation for the time delay between the stress change attending an earthquake and the resultant aftershocks, and even explains why the rate of aftershocks decays as 1/time. Dieterich (31) investigated the behavior of a population of patches having a rate- and state-dependent friction rheology when they are subjected to loading with a periodic component. He found that the response to tidal loading should not be detectable in a homogeneous crust at normal stresses exceeding 8 MPa (or equivalently, at depths greater than about 300 m). It is not clear that any enhancement in correlations associated with the approach to criticality will not likewise be obscured by delayed nucleation. It should be noted that the “lattice solid model” does not yet include rate- and state-dependent friction and the associated delayed nucleation, and thus the results of Wang *et al.* (28) can not be taken as theoretical verification of the LURR effect.

Ouillon and Sornette (32) have performed tests of the concept that large earthquakes are “critical points” on rockbursts in deep South African mines. First, using the simplest signature of criticality in terms of a power law time-to-failure formula, they find evidence both for accelerated seismicity and for the presence of log-periodic behavior [initially proposed by Sornette and Sammis (15)] in the cumulative Benioff strain with a preferred scaling factor close to 2. They have also proposed a new algorithm based on a space and time smoothing procedure, which is also intended to account for the finite range and finite duration of mechanical interactions between events. This new algorithm provides a much more robust and efficient construction of the optimal correlation region, which allows them the use of the log-periodic formula directly in the search process. A preliminary test using the largest event in the catalog shows a remarkable improvement in accuracy and robustness of the perdiction. This suggests new potential for improving on the simple power law fits by using log-periodic signals.

Until recently, this critical point hypothesis has been a conceptual model based on the analogy with phase transitions. Theoretical support has come from simple models such as cellular automata, with (refs. 33 and 34; M. Anghel, W. Klein, J. B. Rundle, and J. S. Sá Martins, unpublished data) and without (13, 14) long-range interaction and granular simulators (ref. 35; P. Mora and D. Place, unpublished data). Models of regional seismicity with more faithful fault geometries have been developed that also show accelerating seismicity before large model events [behavior compatible with the critical point hypothesis (refs. 36–41; G. C. P. King and D. D. Bowman, unpublished data)].

In their cellular automaton model of earthquake faults, Anghel *et al.* (M. Anghel, W. Klein, J. B. Rundle, and J. S. Sá Martins, unpublished data) find that the scaling of earthquake events in models of faults with long-range stress transfer is composed of at least three distinct regions, corresponding to three classes of earthquakes with different underlying physical mechanisms. The largest “breakout” events are found to be associated with a spinodal critical point. In this model, Sá Martins *et al.* (34) find a smoothing effect of the dynamics on the stress field in the precursory phase before a large event.

Huang *et al.* (13) and Sammis and Smith (14) have studied a simple cellular automaton model of earthquakes on a preexisting hierarchical fault structure. The system is found to self-organize at large times in a stationary state with a power law Gutenberg–Richter distribution of earthquake sizes. The largest fault carries irregular great earthquakes preceded by precursors developing over long time scales and followed by aftershocks obeying an Omori's law. The cumulative energy released by precursors follows a time-to-failure power law with log-periodic structures, qualifying a large event as an effective dynamical (depinning) critical point. Down the hierarchy, smaller earthquakes exhibit the same phenomenology, albeit with increasing irregularities.

Mora *et al.* (35) and Mora and Place (unpublished data) have used the lattice solid model (42) to describe discontinuous elasto-dynamic systems subjected to shear and compression. They find accelerating seismic energy release in the lead-up to large earthquake events. They also document strong evidence that the stress–stress correlation length grows stronger as a large model earthquake is approached, in agreement with the prediction of the critical point concept (P. Mora and D. Place, unpublished data).

Fisher *et al.* (39) and Dahmen *et al.* (40) showed analytically and numerically that the model of Ben-Zion and Rice (36) of a two-dimensional discrete brittle fault system in a three-dimensional elastic solid has an underlying critical point of a second-order phase transition. The critical point in this model is associated with specific values of tuning parameters (dynamic weakening and conservation of stress transfer during failure events) so the model behavior is not SOC.

Ben-Zion (37) used a modified version of the model of Ben-Zion and Rice (36) that includes both brittle and creep processes, and showed that large model earthquakes are associated with nonrepeating cyclical establishment and destruction of long-range stress correlations on the fault. During gradual tectonic loading, small and intermediate earthquakes produce stress roughening over their size scales, which collectively smooth the longer wavelength components of stress and develop long-range stress correlations on the fault. The pattern is reversed during large ruptures of size approaching the system dimension, which reduce the stress level and smooth the fluctuations along the large rupture area, while creating large stress concentrations near its boundary and increasing considerably the stress outside it. These system-size events reroughen the long wavelength stress field on the fault, destroy the long-range correlation in the system, and set the beginning of a new large earthquake cycle.

In a model of heterogeneous faults interacting through elastic coupling, Heimpel (38) has found the existence of characteristic (but irregular) quake cycles that are broken up into four stages in agreement with the critical point hypothesis: (*i*) relaxation with aftershocks following the main-shock sequence of the previous cycle; (*ii*) self-organization in which stress and strength increase; (*iii*) criticality characterized by large stress/strength fluctuations; and (*iv*) a main shock that is accompanied by a rapid stress and strength drop.

Using a regional lithospheric model consisting of a seismogenic upper crust governed by damage rheology over a viscoelastic substrate, Ben-Zion and Lyakhovsky (41) demonstrate the existence of accelerated seismic release prior to large model earthquakes, in the case when the seismicity preceding a given large earthquake has a broad frequency-size statistics. The accelerated seismicity is found to be accommodated both by increasing rates of moderate events and increasing average event size.

G. C. P. King and D. D. Bowman (unpublished data) developed a model that is based solely on the decay of the stress shadow from a previous large event. In their model, the stress shadow, calculated using an elastic dislocation model of a great earthquake, is relaxed linearly in time. At each time step, fractal noise is added to the stress and events are calculated for those areas above a failure threshold. Increases in event size and accelerated seismic release are produced by an increase in the number and size of patches above the failure threshold as the shadow decays.

Although each of these models appears to simulate observed seismic acceleration before great events, their formulations are quite different and their relation to the critical point hypothesis is not always clear. The works of Fisher *et al.* (39) and Dahmen *et al.* (40) develop clear relations between seismicity and criticality, although they do not address the issue of accelerating seismicity. In this paper, we attempt to clarify the critical point concept by pointing out that a critical point occurring as a function of time can be better seen as a so-called “finite-time singularity,” a mathematical concept encompassing and generalizing that of a “critical point.” Another advantage of this terminology is that “critical point” is loaded with many different meanings and has been used in distinct ways in several contexts and scientific communities. We exemplify the finite-time singularity concept by simple mechanical models and then explore how models proposed to simulate regional seismicity give rise to a finite-time singularity.

## Finite-Time Singularities

The mathematics of singularities is applied routinely in the physics of phase transitions to describe, for instance, the transformations from ice to water or from a magnetized to a demagnetized state when raising the temperature, as well as in many other condensed matter systems. Such singularities characterize so-called critical phenomena and critical points. This was the concept proposed initially by Sornette and Sammis (15) for the description of large earthquakes. In these problems, physical observables such as susceptibilities, specific heat, etc., exhibit a singularity as the control parameter (temperature, magnetic field) approaches a critical value.

It may be useful to be more precise and view large earthquake as belonging to a larger class of singularities occurring in dynamical systems, which are spontaneously reached in finite time. Because earthquakes are organized in time with a great event terminating a “cycle,” this seems to be a useful view point. Spontaneous singularities are quite common and have been found in many well established models of natural systems, either at special points in space such as in the Euler equations of inviscid fluids (43, 44), in the surface curvature on the free surface of a conducting fluid in an electric field (45), in vortex collapse of systems of point vortices, in the equations of General Relativity coupled to a mass field leading to the formation of black holes (46, 47), in models of microorganisms aggregating to form fruiting bodies (48), or in the more prosaic rotating coin (Euler's disk; ref. 49).

The simplest representative dynamical evolution equation leading to a finite-time singularity is Its solution is where the critical time *t*_{c}*= (m − 1)/[E(0)] ^{m−1}* is determined by the initial condition

*E(0).*The singularity results from the fact that the instantaneous growth rate

*d*ln

*E/dt = E*is increasing with E and thus with time. This can be visualized by studying the “instantaneous” doubling time, defined at the time interval

^{m−1}*Δt*necessary for

*E(t)*to double—i.e.,

*E(t + Δt) = 2E(t),*if the growth rate were fixed constant at its value at time

*t.*When the growth rate of

*E*increases as a power law of

*E,*the doubling time decreases fast and the sequence of doubling time intervals shrinks to zero sufficiently fast so that its sum is a convergent geometrical series. The variable thus undergoes an infinite number of doubling operations in a finite time, which is the essence of a finite-time singularity. The growth of the growth rate captures a positive feedback of the variable

*E*on its growth rate. This positive feedback is the main ingredient for a “finite-time singularity.” In the context of earthquake nucleation, Dieterich (ref. 16; see ref. 8 for a simple derivation) has shown that the usual slip- and velocity-dependent solid friction law introduced in the equation of elasticity gives rises to a finite-time singularity of the slip velocity, representing the transition from accelerated nucleating slip and the fast elasto-dynamical rupture regime. In this case, the positive feedback results from the fact that as the slip distance and slip velocity increase, the solid friction coefficient decreases, thus enhancing further slip.

We now discuss several other examples in which positive feedback leads to a finite-time singularity. We begin with models of crack growth and then move to models with distributed failure that are more suitable for regional seismicity.

## Positive Feedback in Crack Growth

A well known example of a positive feedback occurs in the subcritical crack growth usually represented by a power law equation where *L* is the crack length, *K* is the stress intensity factor, and *C* and *p* are material parameters. In slow early deformation phases, the stress corrosion exponent is typically in the range 2–5. Using the standard result from continuum elasticity *K ∝ L ^{1/2},* the subcritical crack growth equation becomes Because

*p/2 > 1,*this is nothing but Eq. 2 with

*m = p/2,*which gives the finite-time singularity for the crack length according to Eq. 3. This singular behavior has been noted by several authors and has been used for instance to explain the power law distribution of fault length (50). Ben-Zion and Lyakhovsky (41) then convert it into an equation of the form of Eq. 1 by using the relationship between seismic moment and fault length

*M*

_{0}∝ L^{2−3}.The finite-time singularity of crack growth has recently been refined (S. Gluzman and D.S., unpublished data) within a self-consistent theory of crack growth controlled by a cumulative damage variable *d(t),* dependent on stress history in an elastic medium, in the quasistatic regime where the sound wave velocity is taken as infinite. Depending on the damage exponent *m,* which controls the rate of damage *dd/dt ∝ σ ^{m}* as a function of local stress σ, two regimes are found. For

*0 < m < 2,*the model predicts a finite-time singularity. For

*m ≥ 2,*the rupture dynamics require a regularization scheme in the form of a saturation of damage, a minimum distance of approach to the crack tip, or a fixed stress maximum. This shows that the theory leading to a strong finite-time singularity has no continuous limit and is controlled by microscopic processes close to the crack tip. We stress that the mathematical existence of a finite-time singularity is nothing but the signature of a transition to another regime—here the transition to an elasto-dynamic rupture controlled by the finiteness of the speed of elastic shear waves.

Gluzman *et al.* (51) have used insight on the existence of finite-time singularities and their associated power law behavior to develop new theoretical formulas for the prediction of the rupture of systems that are known to exhibit a critical behavior, based solely on the knowledge of the early time evolution of an observable, such as the precursory acoustic emission rate or seismic rate as a function of time or of stress. From the parameterization of such early time evolution in terms of a low-order polynomial, the functional renormalization approach transforms this polynomial into a function that is asymptotically a power law. The value of the critical time *t*_{c}, conditioned on the assumption that it exists, can be determined from the knowledge of the coefficients of the polynomials. This prediction scheme shows a much stronger robustness and reliability with respect to the order of the polynomials and as a function of noise compared to direct power law fits.

## Positive Geometrical Feedback in Laboratory Creep Rupture

Consider the problem of so-called creep rupture (52) in which a rod is subjected to uniaxial tension by a constant applied axial force *P.* The cross section *A(t)* of the rod is assumed to be a function of time. The problem is simplified by assuming that *A(t)* is independent of the axial coordinate, which eliminates necking as a possible mode of failure. The creep deformation is assumed to be isochoric—i.e., the rod volume remains constant during the process. This provides a geometric relation between the cross-sectional area and length *A _{0}L_{0} = A(t)L(t) =* constant, which holds for all times.

The rate of creep strain *ɛ _{c}* can be defined as a function of geometry as showing that

*ɛ*ln

_{c}=*(L(t)/L*ln

_{0}) = −*(A(t)/A*where

_{0}),*L*and

_{0}= L(t_{0})*A*correspond to the undeformed state

_{0}= A(t_{0})*ɛ*at time

_{c}= 0*t*Assuming power-law creep where

_{0}.*σ = P/A*is the ratio of the applied force over the cross-section of the rod. Eliminating

*dɛ*between Eqs. 6 and 7 and using

_{c}/dt*σ = P/A*leads to

*A*the solution of which is

^{μ−1}(dA/dt) = −CP^{μ},*A(t) = A(0)[(t*

_{c}

*− t)/t*

_{c}

*]*where the critical failure time is given by

^{1/μ},*t*

_{c}

*= [A(0)/P]*The rod cross section thus vanishes in a finite time

^{μ}/(μC).*t*

_{c}

*− t*and as a consequence the stress diverges as the time

_{0}*t*goes to the critical time

*t*

_{c}as

*σ = P/A = [P/A(0)][(t*

_{c}

*− t)/t*

_{c}

*]*Physically, the constant force is applied to a thinner cross section, thus enhancing the stress, which in turn accelerates the creep strain rate, which translates into an acceleration of the decrease of the rod cross-section and so on. In other words, the finite-time singularity results from the positive feedback of the increasing stress on the thinner geometrical cross-section and vice versa. This finite-time singularity for the stress can be reformulated as a self-contained equation expressed only in terms of the stress:

^{−1/μ}.*dσ/dt = Cσ*which is indeed of the form

^{μ+1},**1**because μ > 0.

## Positive Feedback in a Simple Damage Model

Sammis *et al.* (53) demonstrated that Eq. 1 follows from the most elementary form of damage mechanics. Their argument is closely analogous to that in the previous section and goes as follows. Let *V _{i}(t)* be the volume fraction of intact rock at time

*t*and

*V*

_{d}

*(t)*be the volume fraction of damage rock such that

*V*

_{i}

*(t) + V*

_{d}

*(t) = 1.*If the damaged rock is not able to support load, then the average local stress, σ

_{local}, will be related to the regional remote stress, σ

_{remote}, as

*σ*

_{local}

*= (σ*

_{remote}

*/A*

_{i}

*) = (σ*

_{remote}

*/V*

_{i}

*)*(because the average cross-sectional area fraction equals the volume fraction). If the average rate of damage production is proportional to the average local stress raised to some power β > 0 (Norton's law), we have Eq. 8 can be integrated to yield

*t*

_{c}

*− t = cV*(when

*t ≈ t*

_{c}), where

*V*

_{i}

*(t*

_{c}) = 0. If the cumulative energy released at time

*t, E(t),*is proportional to the damaged volume

*V*

_{d}

*(t)*we have which has the form of Eq. 1. Note that the structure of Eq. 9 remains robust close to

*t*

_{c}if

*E(t)*is an arbitrary monotonously increasing function of

*V*

_{d}

*(t).*The finite-time singularity at

*t = t*

_{c}in Eq. 9 is the result of positive feedback in Eq. 8: because

*σ*

_{local}

*= (σ*

_{remote}

*/V*

_{i}), we get

*dσ*

_{local}

*/dt ∝ σ*, which is indeed of the form of Eq. 2 with β > 0. The positive feedback comes from the geometrical feedback of the damage onto the local stress.

## Positive Feedback in a Damage Model for Regional Seismicity

Ben-Zion and Lyakhovsky (41) calculate regional seismicity by using continuum damage rheology. In their model damage creation is linked to the strain. In one dimension it takes the form *dα/dt = cɛ ^{2},* where α is the damage and ɛ is the one-dimensional strain. Feedback in this case comes through the elastic constant

*E = E*which is an explicit function of the damage. Hooke's law is thus

_{0}(1 − α)*σ = E*which, at constant stress, leads to In Eq. 10, α = 1 when

_{0}(1 − α)ɛ,*t = t*

_{c}which defines the critical time. Although the strain has a singularity at

*t = t*

_{c}, the cumulative Benioff strain does not. Following Ben-Zion and Lyakhovsky (41), the elastic energy release rate can be written which, when integrated to give the Benioff strain, gives consistent with the observed values. This model is another example of a finite-time singularity caused by positive feedback. Differentiating Hooke's law and using

*dα/dt = cɛ*gives

^{2}*(dɛ/dt) = (cE*which has the form of Eq. 2 and leads to the finite-time singularity in ɛ.

_{0}/σ)ɛ^{4},## Positive Feedback in a Percolation Model of Regional Seismicity

In the previous examples, the finite-time singularity was the result of positive feedback in the failure process at constant stress. It is also possible to have a finite-time singularity in a heterogeneous system driven by increasing stress, as when the stress shadow from a prior great earthquake decays because of tectonic loading. Imagine that the heterogeneity of the crust is accounted for by stochastic rupture thresholds σ_{th}*(r),* where *r* is position, distributed according to some distribution with only short-range spatial correlation. Let us assume that the last greatest earthquake has cast a stress shadow −σ_{0} over a very large area *S* and let us neglect for the time being its spatial dependence and assume that it is uniform in space. Tectonic loading then progressively increases the stress everywhere uniformly at the rate *dσ*_{tect}*/dt.* The stress-to-rupture at an arbitrary point *r* at time *t* is thus We assume that only those points *r* such that *Δσ(r) > 0* are activated and can rupture as earthquakes. The other points with negative stress-to-rupture are inactive. The field *Δσ(r)* can be represented as a highly corrugated random mountainous landscape in which the altitude *Δσ(r) = 0* defines the boundary between seismic active regions *[Δσ(r) > 0]* and inactive domains *[Δσ(r) < 0].* The altitude *Δσ(r) = 0* can be envisioned as a “sea level.” Immediately after the large earthquake, most of the points *r* have *Δσ(r) < 0,* which means that only a few islands emerge above sea level. Ignoring aftershocks, the seismic activity is very weak, localized, and incoherent. As time goes on and tectonic loading increases, the “mountainous landscape” *Δσ(r)* is progressively lifted up. As a consequence, a larger and larger distribution of islands emerges above sea level, until a point when macroscopic “continents” appear. This process is nothing but the percolation model governed by the control parameter defined as the fraction *p* of territory above sea level. At early time after the great event, *p* is very small and the stress landscape is below percolation. As time increases, *p* increases approximately linearly with time, until it reached a critical value *p*_{c}, called the percolation threshold, corresponding to the appearance of very large clusters of connected points above sea level.

Percolation theory (54) provides a precise prediction on the distribution of cluster sizes *n(s)*—i.e., the distribution of earthquake sizes in this simple model. It is known that, sufficiently close to *p*_{c}, the number of clusters of size *s* per unit area is given by where the typical maximum cluster size diverges as as *p* goes to the critical point *p*_{c}. *s*(p)* is nothing but the correlation length of the system. τ and σ are two critical exponents that depend on the dimension of the system. In two dimensions, τ = 187/91 and σ = 36/91, whereas in high dimensions or in mean-field theory τ = 5/2 and σ = 1/2. Note that, at the critical point *p = p*_{c}, the distribution of earthquake size becomes exactly self-similar according to a pure power law. By definition, *∑*_{s}*sn(s) = p.* The average size of earthquakes (clusters) is then given by where γ is another critical exponent, related to τ and σ by the scaling law: γ = (3 − τ)/σ. The released seismic energy is usually taken as the 3/2 power of the rupture area *s.* This shows that the average seismic energy released per earthquake is expected to increase as the power law *|p − p*_{c}|^{−3γ/2} as the critical point is approached. This simple stress percolation model thus predicts an accelerated increase with time of the typical largest earthquake *s*(p)* and of the average seismic energy *〈s〉 ^{3/2}* according to the law of finite-time singularity. Eq. 16 can indeed be rewritten as is

*(d〈s〉/dt) = 〈s〉*with λ = (γ + 1)/γ > 1. This expression can also be recovered from a mapping of the random point percolation problem on a version of the forest fire model using a cluster analysis (55).

^{λ}There is, however, a problem. For random percolation, the power γ is always greater than 1 {e.g., in two dimensions, γ = [(3 − τ)/σ] ≈ 3} so that not only is average seismic energy in a given time interval (i.e., the energy rate) singular, but so is its integral, the cumulative energy (or cumulative Benioff strain). Such is not the case when the percolating topography is correlated. In particular, it is well known that rupture can be modeled as correlated percolation (56, 57): the correlation strength continuously evolves in a system progressively brought to failure by the action of the stress field that provides an additional channel of communication on top of the sole connectivity defined in standard percolation. A conceptual model of such correlated percolation mediated by the existence of an additional field such as the stress is to consider spins interacting with ferromagnetic interactions (playing the role of the stress field) on a percolating lattice. In this case, the exponents can be calculated and one finds that the exponent σ is changed to 2 rather than of the order of 1/2 or 1/3 and the exponent γ becomes less than 1 (ref. 54, p. 146). In this class of correlated percolation models with an addition field, the probability to be above sea level is a function of past events due to stress redistribution. In this sense, stress redistribution gives positive feedback in the percolation model where *p* is not only growing globally as a function of time but is modified locally due to stress redistribution by prior events. For the cumulative Benioff strain, which is proportional to the integral over time of the average cluster size, this correlation in time with path-dependent history leads to γ < 1 and therefore to a finite-time singularity of the form of Eq. 1.

## Positive Feedback in a Stress-Shadow Model for Regional Seismicity

In the previous model, the accelerating seismicity and finite-time singularity emerge as the direct consequences of two ingredients: (*i*) the power law dependence of the correlation length on the distance to the critical point and (*ii*) the sweeping of the control parameter *p,* which increases monotonically with time (see ref. 58 for a general discussion of sweeping of control parameters in critical point models). In this model, the accelerated seismicity is expressing the growth of correlation of stress within the system quantified by the size of the clusters above the rupture threshold. King and Bowman (unpublished data) modify the uniform stress shadow by calculating the spatial redistribution of stress predicted by elastic dislocation theory, to which they added fractal noise. Their tectonic loading is not uniform, but also consists in decreasing progressively the strength of the elastic dislocation equivalent to the great earthquake. In addition, they allow the noise to be redistributed at each time step. When they increased this regional stress field linearly in time to simulate tectonic loading, and applied a threshold failure criterion, their model seismicity was well fit by Eq. 3. In their model, *t*_{c} is the time at which tectonic loading has eliminated the stress shadow. We now give an analytic solution to this problem with the goal of deriving Eq. 1.

Consider the section shown in Fig. 1 orthogonal to the fault on which the great event last occurred. This event is supposed to cast a stress shadow, which, for simplicity, we take as a negative scalar field σ_{sh}*(r)* that is only a function of the distance *r* from the fault: where −σ_{0} is the stress drop on the fault, *h* is the depth of the fault, and *d* is the dimension of faulting *(d = 2* for fracture through a plate, *d = 3* for a half space). To this stress shadow we add random stress fluctuations to represent strength heterogeneity and stress perturbations associated with prior events. We will assume a constant failure threshold. This is equivalent to assuming random failure thresholds and spatially smooth stress level as in the previous section. The second ingredient of the model is a simulation of tectonic loading by decreasing the strength of the dislocation linearly with time according to *σ _{0}(t) = σ_{0}(0) − dσ*

_{tect}

*/dt.*The total average stress applied at a point

*r*from the fault is thus where the origin of time is taken at the time of the last great earthquake. As the third ingredient of the model, we assume that intermediate sized earthquakes only become possible when the average stress level reaches some threshold −σ* < 0 corresponding to the percolation threshold discussed above. At time

*t,*we can thus define a boundary

*r*(t)*such that all points

*r > r*(t)*are above the percolation threshold and can generate intermediate events, while all points

*r < r*(t)*are still in the shadow below the percolation threshold and can only produce small events (see Fig. 1). The transition distance

*r*(t)*is such that which gives where

*t*

_{c}

*≡ σ*

_{0}(0)/(dσ_{tect}

*/dt),*the time at which the tectonic stress has reloaded the main fault so that it can rupture again. Time

*t**is defined at the time at which

*σ(0, t*) = −σ**and the entire shadow has reached the percolation threshold. Hence,

*r*(t*) = 0*(see Eq. 20). We also may write

*t*

_{c}

*− t* = σ*/(dσ*

_{tect}

*/dt).*Note that

*r*(t)*goes to zero in finite time as

*t → t**with a singular behavior having exponent

*1/d.*Thus, new area reaches the percolation threshold at a rate proportional to

*−(dr*/dt).*

If the total area were constant, the cumulative number of earthquakes would simply grow linearly with time. In the present case, the active area grows with time. We assume that when an elementary domain sees its total applied stress *σ(r)* given by Eq. 16 reach the percolation threshold σ*, a sudden burst of earthquake energy is released as a result of the rapidly growing correlation. After this brief burst, the seismic activity assumes a constant background level. As a consequence, the cumulative energy released by earthquakes, *E(t),* obeys the following differential equation where *a* is the constant background seismicity of the region freed from the large earthquake shadow per unit area, *L* is the length of the fault, and *g* quantifies the intensity of the burst of seismic activity triggered by the percolation. *R* is introduced as a convenient regularization, which has the physical meaning of the total size of the region influenced by the great earthquake. The results are independent of its specific value. We stress that the second term *−gL(dr*/dt)* of Eq. 21 is present as soon as the seismicity at the border of the shadow region is distinct from that outside the shadow. This term can for instance take into account an overshoot of the earthquakes suddenly freed and/or hysteresis of the mechanical deformations.

Neglecting the second term on the r.h.s. and for constant *r*,* Eq. 21 becomes *dE/dt = aL(R − r*),* the solution of which is *E(t) = aL(R − r*)t,* corresponding to the usual linear-in-time cumulative number of earthquakes. Close to *t*,* the integration of the first term, *aL(R − r*(t)),* gives a nonsingular term the leading behavior of which is proportional to time *t.* Close to *t*,* it can be approximated by a constant *E _{0}* because it varies much more slowly and more smoothly than the other term. More interestingly, the second term

*−gL(dr*/dt)*gives rise to the singular behavior of Eq. 20 with an infinite accelerating slope close to

*t*

_{c}. Putting everything together, we get finally that the cumulative energy release is governed by It is interesting that stress on the main fault does not fully recover at

*t = t*,*but requires an additional time increment

*t*

_{c}

*− t* = σ*/(dσ*

_{tect}

*/dt)*that depends on the percolation threshold σ*.

Expression **22** has exactly the form of Eq. 1 with *m = 1/d* and *B = −ghL/(t*_{c}* − t*).* This calculation provides possibly the simplest example of how an accelerating power law behavior culminating at a finite-time singularity may emerge from a simple geometrical mechanism. The exponent *m = 1/d* of the singular behavior is solely controlled by the *1/r ^{d}* dependence at large

*r*of the stress shadow σ

_{sh}given by Eq. 14. If

*d = 2*for a large rupture that spans the brittle crust then

*m = 1/2,*if

*d = 3*for a smaller rupture contained in the crust, then

*m = 1/3.*

The functional form of Eq. 22 is unchanged when including explicitly heterogeneous residual stress field and stress thresholds that vary from point to point in a random way. A statistical theory extending the previous calculation can be formulated that amounts to summing over the probability distribution of threshold and over space. The exponent *1/d* is modified into *b + 1/d* if the distribution of stress thresholds is of the form *P(σ*_{th}*) ∝ σ**,* where *0 ≤ b* for small σ_{th}—i.e., if there are few elements with weak rupture thresholds close to 0. The larger *b* is, the stronger is the rock system, because there are fewer and fewer weak elements. Therefore, in this theory, the finite-time singularity relies on the existence of a significant fraction of weak faults or cracks (such that *b* is close to 0) that are destabilized as the stress shadow disappears.

## Discussion

Models that have been proposed to describe accelerating seismicity before large earthquakes fall into two general classes: those based on mechanical feedback and those based on the decay of the stress shadow from a prior event. In the former, the power law acceleration of regional seismicity is due to positive feedback in the failure process at constant large-scale stress. This feedback can be the result of stress transfer from damaged to intact regions, or it can result from the effect of damage in lowering the local elastic stiffness. Both result in finite-time singularities.

Models based on the recovery of a stress shadow due to tectonic loading and stress transfer are quite different. Accelerating seismicity results from the increasing stress as the shadow is eroded by tectonic loading. In this case, a finite-time singularity can result from either correlated percolation (describing stress transfer between intermediate size events) in a uniformly rising stress field, or from the geometrical shrinking of the *r ^{−d}* stress perturbation coupled with random percolation. In the latter case, the finite-time singularity results in a “front” of intermediate sized events that propagates toward the main fault as the area closer to the fault progressively reaches the percolation threshold.

Although both classes of models predict accelerating seismicity and growth of event size before a great earthquake, they predict different spatial patterns for this precursory seismicity that may offer an observational way to distinguish between them. Models based on feedback in the failure process at constant stress predict rapid shear localization with the largest precursory event occurring in or near the fault plane of the ultimate event. Models based on the decay of stress shadows, on the other hand, predict that the largest precursors should develop first at the fringes of the shadow and gradually move toward the plane of the great event. Current observations tend to find intermediate precursory events at large distances (22, 24), but a quantitative test of the two classes of models has yet to be performed.

Finally, we end with a somewhat speculative figure that tends to support a mechanism based on the decay of stress shadows. Amos Nur (personal communication) has observed that strike slip earthquakes in the Mojave Desert have been occurring progressively closer to the San Andreas Fault since 1857. In Fig. 2, the ordinate is the distance of each event from the closest point on the San Andreas Fault, and the abscissa is the elapsed time since the last great event on that portion of the fault (1857 to the north, 1812 to the south). Coordinates of the San Andreas fault were taken from the U.S. Geological Survey data base at http://geohazards.cr.usgs.gov/eq/html/faults.html. Earthquake epicenters and magnitudes were taken from the *Ellsworth Catalog* of historical seismicity that can be found at http://pasadena.wr.usgs.gov/info/cahist_eqs.html, and from the map at the Southern California Earthquake Center (SCEC) data center (http://www.scecdc.scec.org/clickmap.html). Although there are not many events, there is a monotonic decrease in distance with time.

## Acknowledgments

This work was supported by National Science Foundation Grants EAR9804870 (to C.G.S.) and DMR9971475 (to D.S.), and by the James S. McDonnell Foundation 21st century scientist award/studying complex system (to D.S.).

## Footnotes

↵† To whom reprint requests should be addressed. E-mail: sammis{at}usc.edu.

This paper results from the Arthur M. Sackler Colloquium of the National Academy of Sciences, “Self-Organized Complexity in the Physical, Biological, and Social Sciences,” held March 23–24, 2001, at the Arnold and Mabel Beckman Center of the National Academies of Science and Engineering in Irvine, CA.

## Abbreviations

SOC, self-organized criticality

LURR, Load/Unload Response Ratio

- Copyright © 2002, The National Academy of Sciences

## References

- ↵
- Sieh K.

**,**603-623. - ↵
- Bakun W. H.

**,**3051-3058. - ↵
- Nur A.

**,**1217-22. - Whitcomb J. H.

**,**632-641.- Scholz C. H.

**,**803-809.- Griggs D. T.

**,**537-540.- ↵
- ↵
- Sornette D.

- ↵
- ↵
- Geller R. J.

**,**1616-1617. - ↵
- ↵
- ↵
- Huang Y.

**,**43-48. - ↵
- ↵
- Sornette D.

**,**607-619. - ↵
- ↵
- Harris R. A.

**,**229-232. - ↵
- ↵
- Sykes L. R.

**,**595-599. - ↵
- ↵
- ↵
- ↵
- ↵
- Jaumé S. C.

**,**279-306. - ↵
- ↵
- ↵Yin, X-C. (2001)
*Pure Appl. Geophys.*, in press. - ↵
- Wang Y.

- ↵
- ↵
- ↵
- ↵
- ↵
- Anghel M.

**,**F923-F924. - ↵
- ↵
- Mora P.

- ↵
- Ben-Zion Y.

**,**14109-14131. - ↵
- Ben-Zion Y.

**,**5677-5706. - ↵
- Heimpel M.

**,**865-868. - ↵
- ↵
- ↵Ben-Zion, Y. & Lyakhovsky, V. (2001)
*PAGEOPH*, in press. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Choptuik M. W.

**,**353-365. - ↵
- ↵
- ↵
- Sornette D.

**,**1079-1081. - ↵
- Gluzman S.

**,**122-137. - ↵
- Krajcinovic D.

- ↵
- Sammis C. G.

- ↵
- Stauffer D.

- ↵
- ↵
- Sahimi M.

**,**214-395. - ↵
- Johansen A.

**,**163-181. - ↵
- Sornette D.

**,**209-221.

## Citation Manager Formats

## Sign up for Article Alerts

## Jump to section

- Article
- Abstract
- Finite-Time Singularities
- Positive Feedback in Crack Growth
- Positive Geometrical Feedback in Laboratory Creep Rupture
- Positive Feedback in a Simple Damage Model
- Positive Feedback in a Damage Model for Regional Seismicity
- Positive Feedback in a Percolation Model of Regional Seismicity
- Positive Feedback in a Stress-Shadow Model for Regional Seismicity
- Discussion
- Acknowledgments
- Footnotes
- Abbreviations
- References

- Figures & SI
- Info & Metrics