Under the Hood of the Earthquake Machine: Toward Predictive Modeling of the Seismic Cycle

See allHide authors and affiliations

Science  11 May 2012:
Vol. 336, Issue 6082, pp. 707-710
DOI: 10.1126/science.1218796


Advances in observational, laboratory, and modeling techniques open the way to the development of physical models of the seismic cycle with potentially predictive power. To explore that possibility, we developed an integrative and fully dynamic model of the Parkfield segment of the San Andreas Fault. The model succeeds in reproducing a realistic earthquake sequence of irregular moment magnitude (Mw) 6.0 main shocks—including events similar to the ones in 1966 and 2004—and provides an excellent match for the detailed interseismic, coseismic, and postseismic observations collected along this fault during the most recent earthquake cycle. Such calibrated physical models provide new ways to assess seismic hazards and forecast seismicity response to perturbations of natural or anthropogenic origins.

Seismic and geodetic observations provide an increasingly detailed insight into fault motion over a wide range of temporal and spatial scales, from rapid seismic rupture to slower postseismic slip and complex interseismic behavior, including slow episodes of accelerating slip and tremor (15). Laboratory experiments and theoretical developments (613) provide an increasingly detailed physical basis for understanding the entire earthquake cycle. Yet, models capable of capturing a wide range of observations are still in their infancy. Existing models are either restricted to specific aspects of fault behavior (e.g., progression of a single dynamic rupture or evolution of postseismic slip) or simplify some stages of the fault deformation (1420). Recently developed numerical methods (2123) allow us to resolve, in one physical model, slow tectonic loading, earthquake nucleation, and rupture propagation—including the radiation of seismic waves—and the afterslip transient that follows main shocks, but so far these methods have been applied to qualitative studies of conceptual fault scenarios (24).


Video Lab

We have developed a fully dynamic model of the earthquake cycle capable of quantitatively reproducing a wide range of observations for the Parkfield segment of the San Andreas Fault (SAF) (Fig. 1). The Parkfield sequence of moment magnitude (Mw) 6.0 events, their inferred similarities, and their short recurrence times inspired one of the most famous prediction experiments and prompted the installation of modern seismic and geodetic networks (25, 26). The latest rupture of 2004 defied the expectations by taking place a decade later than anticipated and initiating on the opposite end of the segment compared with previous events (27). Interestingly, a series of smaller earthquakes occurred in 1993 around the location where the 1966 event had nucleated (Fig. 1), but they failed to generate the Mw 6.0 event that was expected at the time (1). Due to the dense instrumentation networks and other observational facilities, such as the San Andreas Fault Observatory at Depth (13), installed to monitor the Parkfield segment, the pattern of strain buildup, microseismic activity in the interseismic period, as well as co- and postseismic deformation related to the 2004 earthquake cycle, have been relatively well constrained (2731).

Fig. 1

(A) Tectonic setting and seismic cycle at Parkfield. The SAF accommodates most of the relative motion between the Pacific (PAC) and North-American (N-AM) plates in central California and produces localized microseismicity (black dots). Deformation during the co-, post-, and interseismic periods is monitored by various arrays of instruments, including GPS (black triangles and pre-2004 velocity vectors) and broad-band seismometers (gray squares). (B) Spatial distribution of fault slip and microseismicity on the SAF during the earthquake cycle. Aseismic creep occurs above and below the SZ. The area between the microseismic streaks is coupled. The coseismic slip distribution of the 2004 main shock (20-cm red contours) may be concentrated in the domain circumscribed by background seismicity. (C) Paleoseismic cycle of the largest earthquakes since 1881 (with interevent times between 12 and 38 years). Smaller events occurred in 1992 to 1994, culminating with the November 1993 Mw 4.6 and the December 1994 Mw 4.7 earthquakes.

Our dynamic model of the earthquake cycles at Parkfield is constrained by multiple sets of observations and previous theoretical findings. We use the spatial pattern of microseismicity, the time series of Global Positioning System (GPS) displacements in the 1999 to 2010 period, the interferometric synthetic aperture radar (InSAR) data, and the GPS offsets of the 2004 earthquake (2, 4, 19, 29, 32). We also consider the slip distribution of the 1966 event and the historical catalog of recurrence times and hypocenter locations of Mw 6.0 events (25). As an integration device for all observations, we adopt a strike-slip fault segment embedded into an elastic medium, loaded by a deep-seated slip at the long-term slip rate, and governed by rate-and-state friction, a well-established empirical constitutive law for fault strength (68, 10). The area with rate-weakening friction, where seismic slip can nucleate, is bounded to the north and south by rate-strengthening patches. The northern rate-strengthening patch accounts for the creeping segment of the SAF. The southern one is more speculative; it serves as a possible proxy for the kind of barrier effect needed to account for the repetition of similar events arresting in that area and for an as-yet-unknown source of localized stressing (35).

The model results in a rich history of fault slip with spontaneous nucleation and ruptures of earthquakes of magnitudes ranging from Mw 2.0 to 6.0 (Fig. 2 and movie S1). The simulated Mw 6.0 earthquake cycles reproduce co-, post-, and interseismic behavior of the Parkfield segment, with most coseismic slip occurring in the area circumscribed by microseismicity. The simulated sequence of earthquakes includes the nucleation of a rupture near the hypocenter of the 2004 Mw 6.0 event (Fig. 2B). The rupture propagates northward, slows down in the vicinity of the velocity-strengthening obstacle (Fig. 2C), and proceeds to rupture the entire seismogenic zone (SZ), arresting at the creeping segment (Fig. 2D), similarly to what occurred during the 2004 event. The resulting seismic moment corresponds to a Mw 6.0 earthquake. During the rupture, the shallow asperities slip with aseismic sliding velocity below 0.1 m/s. After the coseismic rupture, the SZ becomes locked, and afterslip expands around it (Fig. 2E). The afterslip slows down with time, eventually leading to the pattern of the locked SZ surrounded by interseismic creep (Fig. 2F). During the interseismic period, creep penetrates into the SZ, building conditions for earthquake nucleation. Eventually, the cycle repeats, with a mean recurrence time of ∼20 years.

Fig. 2

Model that reproduces the entire seismic cycle at Parkfield. (A) Spatial distribution of rate-weakening (white) and rate-strengthening (gray) friction properties (ab), with rate-weakening values in the SZ delineated by the background seismicity, a few shallow patches required to explain near-field GPS data, and two deeper patches included for illustration purposes. The critical size h* for unstable slip is shown in black circles (see supplementary text). (B to F) Snapshots of a Mw 6.0 seismic cycle with rupture nucleating spontaneously to the south near the Parkfield 2004 earthquake hypocenter (B), propagating to the north and rupturing the entire SZ [(C) and (D)], and followed by a slow postseismic transient (E), with interseismic loading of the partially locked SZ (F). Another Mw 6.0 event nucleates 20 years later. Zero time is chosen for plotting convenience. The solid gray profiles indicate the contours of the cumulative slip at 0.1-m intervals.

Our model explains the longer-term behavior of the earthquake cycle at Parkfield. Some of the events produced offer an excellent fit to the coseismic and the postseismic displacements measured at the GPS stations, indicating an adequate representation of faulting over a range of time scales (Fig. 3 and fig. S3). It also qualitatively reproduces the switch in the hypocenter location from earthquakes nucleating at the northern end to earthquakes nucleating at the southern end, as occurred in 2004 at Parkfield (Fig. 4). The southern and northern corners of the SZ are favorable nucleation sites because of their location near local stress concentrations due to both horizontal and vertical boundaries between domains of stable and unstable sliding. The domain with stable sliding at the northern end of the SZ is clearly seen in geodetic inversions; no equivalent domain is observed to the south (Fig. 1B). It is possible that there is a local zone of creep not resolvable with the data used in these inversions but sufficient to induce a local reloading, at least in the early postseismic period. In our model, the transition between the two nucleation sites occurs after a few smaller earthquakes (Mw 2.0 to 4.0) that can be interpreted as failed nucleation attempts of the main event. The transition in nucleation location is coincidental with the longest simulated recurrence time (Trmax = 23.1 years, compared with the smallest, Trmin = 15 years) of the sequence (Fig. 4 and movie S2).

Fig. 3

(A) Comparison between the observed surface displacements of the 2004 Mw 6.0 Parkfield earthquake and those predicted by our numerical simulation. (B) Sample time series of postseismic displacements during a 1-year period after the earthquake: GPS time series at station MASW (black dots) and numerical simulations (red solid line). (C) Map view of postseismic transients in a 2-year period at the Southern California Integrated GPS Network stations: GPS data (circles) and model (squares).

Fig. 4

Comparison between modeled and inferred slip history at Parkfield. (A) Depth profiles in the middle of the SZ comparing the modeled history of fault slip for three events to the one inferred from our inversion of InSAR and GPS data in the period 1999 to 2010. The red stars indicate the location of hypocenters. The seismicity concentrates around the area of maximum coseismic slip. (B) Horizontal profiles in the middle of SZ, at 7.5-km depth, showing a series of 11 simulated ruptures. The series captures a complex transition of hypocenter location from one side of the SZ to another, punctuated by foreshocks and a seismic crisis consisting of Mw 2.0 to 4.0 earthquakes and giving rise to the longest recurrence time of a Mw 6.0 quake. The red dashed profiles represent cumulative coseismic slip at 1-s intervals; the blue solid profiles are cumulative aseismic slip at 1-year intervals taking place during the post- and interseismic periods.

In view of these results, the delay of the latest Mw 6.0 2004 event at Parkfield and its unexpected change of hypocenter may have been at least partially caused by the smaller earthquakes in 1992 to 1994, a series of arrested Mw 6.0 nucleations. On 20 October 1992, close to the expected time of a Mw 6.0 earthquake, a Mw 4.6 event occurred near the 1966 hypocenter (33) as a part of an accelerated slip episode between 1992 and 1994 (1). The Mw 4.6 event failed to grow to a Mw 6.0 event, indicating that the surrounding area was not yet ready to sustain rupture propagation due to cycle-to-cycle variations inherent in such a nonlinear system, outside perturbations, or finer-scale heterogeneity. The stress drop associated with the Mw 4.6 event may have prevented nucleation of the next Mw 6.0 event in that area.

Our results indicate that the Parkfield anomalies in the earthquake cycle can at least partially result from a spontaneous behavior of rate-and-state faults after a complex sequence of smaller interseismic events. The model succeeds in producing returning Mw 6.0 events irregularly spaced in time, which we attribute to stress redistributions induced by smaller events. However, the distribution of the recurrence time in our simulation, between 15 and 25 years, is not as wide as the observed variability at Parkfield (fig. S4). Reproducing the full spread of recurrence times since 1881 may require taking into account additional aspects of fault dynamics, including stress transfers from remote earthquakes and their postseismic transient, the effect of deeper crustal processes (5, 34), or more sophisticated friction laws (11, 12). Representing coseismic fault behavior by the standard rate-and-state formulation with rate weakening is an effective but simplified approach, akin to the typical practice of using linear slip weakening laws to simulate dynamic rupture (15, 17).

The presented study demonstrates the possibility of creating comprehensive physical models of fault zones that integrate geodetic and seismological observations for all stages of the earthquake cycle. As computational resources and methods improve, more realistic fully dynamic simulations will be able to incorporate near-fault damage, local fault nonplanarity, and additional coseismic weakening mechanisms. Such simulations could in principle be used to assess the full range of earthquake patterns that a particular fault system might produce or to assimilate observation about past earthquakes and interseismic loading to assess future seismicity (20). The resulting region-specific models may then reveal how the entire earthquake cycle is affected by various additional factors, such as static and/or dynamic triggering and human intervention in the form of geothermal energy harvesting and CO2 sequestration.

Supplementary Materials

Materials and Methods

Supplementary Text

Figs. S1 to S4

References (3674)

Movies S1 and S2

Data File S1

References and Notes

  1. Acknowledgments: This study was supported partly by the National Science Foundation (grant EAR 0548277 to N.L.), the Gordon and Betty Moore Foundation, and the Southern California Earthquake Center (SCEC). This is Tectonics Observatory contribution #194 and SCEC #1600. Numerical simulations for this study were carried out on the CITerra Dell cluster at the Division of Geological and Planetary Sciences of the California Institute of Technology. J.P.A. worked on this study while on sabbatical at GFZ (German Research Centre for Geosciences) (Potsdam, Germany) with support from the Humboldt Foundation. GPS data used in this study are available from the Scripps Orbit and Permanent Array Center. Numerical data are available upon request to the authors.

Stay Connected to Science

Navigate This Article