Back to the deformation and Stress Change Modeling home page

Back to the team's online recent papers pages

Back to Ross Stein's web page

S. Toda and R. S. Stein,

Response of the San Andreas Fault to the 1983 Coalinga-Nuñez Earthquakes: An Application of Interaction-based Probabilities for Parkfield, J. Geohys. Res., submitted 19 January 2001

 

Response of the San Andreas Fault to the
1983 Coalinga-Nuñez Earthquakes:
An Application of Interaction-based Probabilities for Parkfield

Shinji Toda

Earthquake Research Institute, University of Tokyo, Japan

Ross S. Stein

U.S. Geological Survey, Menlo Park, California

Abstract. The Parkfield-Cholame section of the San Andreas fault, site of six M~6 shocks since 1857 and an unfulfilled earthquake forecast in 1985 [Bakun and Lindh, 1985], is the best monitored section of the world’s most closely watched fault [Roeloffs and Langbein, 1994]. In 1983, two large earthquakes, the M=6.5 Coalinga and M=6.0 Nuñez events, struck 26 km northeast of Parkfield. Seismicity rates climbed along the creeping section of the San Andreas north of Parkfield, and dropped along the locked section to the south. Right-lateral creep slowed or reversed from Parkfield south. Here we calculate that the Coalinga sequence increased the shear and Coulomb stress on the creeping section, causing the rate of small shocks to rise for about 18 months until the added stress was shed by additional fault slip. But the 1983 events decreased the shear and Coulomb stress on the locked Parkfield segment, causing surface creep and seismicity rates to drop for more than 6 years. Creep and aftershock observations suggest that the San Andreas takes up to ten times longer to recover from a stress perturbation in the locked section than where it creeps. We use these observations to cast the likelihood of a Parkfield earthquake in terms of an interaction-based probability, which includes both the renewal of stress following the 1966 Parkfield earthquake and the stress transfer to the San Andreas from the Coalinga events. We estimate that the 1983 shocks reduced the 10-year probability of a M~6 Parkfield earthquake by 22% (from 54±22 % to 42±23%), and that the probability did not recover to its pre-1983 level until about 1991. Thus perhaps we can rationalize why the Parkfield earthquake did not strike in the 1980’s, but not why it was absent in the 1990’s. During the next decade (2001-2011), we calculate a 58±17% probability of a Parkfield earthquake.

1. Introduction

Recent efforts to explain features of earthquake occurrence and probability by stress transfer are predicated on an association between Coulomb stress changes and seismicity rate changes. But while a correspondence between calculated stress increases and observed seismicity rate increases has been reported (e.g., [Reasenberg and Simpson, 1992; Stein, 1999; Toda et al., 1998]), it is difficult to measure seismicity rate decreases on all but faults with the highest rates of microearthquakes. In addition, probability models incorporating earthquake interaction [Dieterich and Kilgore, 1996; Parsons et al., 2000] suffer from the need to average stress changes and constitutive parameters over fault surfaces much larger than the site of earthquake nucleation. Study of the Coalinga-Parkfield interaction helps to overcome these obstacles: First, the 1983 events were large enough to impart significant stress to the San Andreas, but far enough away that the unknown details of the 1983 fault slip have a negligible impact on the stress transfer (Fig. 1). Second, the creeping section has among the highest rates of microearthquakes of any fault in the U.S., so seismicity rate decreases are readily detected [Poley et al., 1987; Wyss et al., 1990]. Third, seismic, geodetic, and surface creep data permit an assessment of the parameters needed for an interaction probability analysis. Finally, the 1934 and 1966 Parkfield shocks nucleated within a 6 x 6 km fault patch [Bakun and McEvilly, 1979; Bakun and McEvilly, 1984], and the earthquake slip in 1934 and 1966 was similar north of Cholame [Segall and Du, 1993], so one can focus on the site of past–and perhaps future–earthquake nucleation and rupture.

Fig. 1 (a) Shear stress transferred by the 2 May 1983 M=6.5 Coalinga (CO) and 11 June-22 July M=6.0 Nuñez (NZ) earthquakes on vertical right-lateral planes parallel to the San Andreas fault at 8 km depth. Source parameters for CO are from Stein & Ekström (1991) for the coseismic period: 150° strike, 15°W dip, 4.7 m reverse slip, 10 km upper depth, 1.5-4.0 km width; for NZ they are based on Eaton [1990] and Rymer et al. [1990]: 178° strike, 65°E dip, 0.22 m right-lateral and 0.65 m reverse slip, 5.4 km length, 2 km upper and 8.3 km width. Excluded are the 25 October 1982 M=4.8 New Idria (NI) shock because of its small moment, and the 4 August 1985 M=6.0 Kettleman Hills (KH) shock [Ekström et al., 1992] because it is farther from the San Andreas; these have negligible impact on the stress. (b) The observed seismicity rate change is superimposed along the San Andreas fault, also at 8 km depth.

2. Observations of Seismicity Rate Change

The seismicity rate change is superimposed on the shear stress change imparted by the Coalinga-Nuñez shocks in Fig. 1b, and these are plotted as a function of depth in Fig. 2. The seismicity rate change for the post-Coalinga 1.5-year period 830502-841101 (yr/mo/dy) (Fig. 2d) and the 5.5-yr period 841102-90501 (Fig. 2e) are calculated relative to the 3-year pre-Coalinga period 80502-830501. Earthquakes from the Northern California Seismic Network (NCSN) were relocated by Felix Waldhauser using a double-difference algorithm [Waldhauser et al., 1999] and the current NCSN Parkfield velocity model. The seismicity rate is measured for M_L > 1.3 shocks starting in 1980, the magnitude of catalog completeness once a shift in magnitude assignments by the USGS (relative to Berkeley magnitudes) was made in 1978-1980, as identified by Stefan Weimer [Weimer and Wyss, 1997]. The number of M_L > 1.3 earthquakes in disks of 5-km radius with centers spaced 1 km apart is counted for pre- and post-Coalinga periods. The rate change is computed and smoothed with a Gaussian filter for every disk in which there are at least 6 pre-Coalinga shocks (n_min>6) [Matthews and Reasenberg, 1988; Reasenberg and Simpson, 1992]. Because more accurately relocated Parkfield seismicity using waveform cross-correlation for post-1990 shocks shows few, if any, earthquakes in the uppermost 2 km [Rubin et al., 1999], we regard rate changes in the upper 2 km in the 1980-1990 data as unreliable, and they are not considered.

Fig. 2 Calculated normal (a), shear (b), and Coulomb (c) stress change and observed seismicity rate change within 5 km of the San Andreas fault (d-f) associated with the 1983 Coalinga-Nuñez shocks, with earthquakes during the indicated post-Coalinga period superimposed in black. The stress changes are similar to those of Simpson et al. [1988]. Note the association between the calculated shear stress change (b) and the observed seismicity rate change during the first 1.5 yr after the Coalinga sequence (d). The rate increase in the creeping section (km = —50 to 0) disappears after 1-2 years (e). As a visual guide to the rate change calculations in (d), seismicity during the 1.5-yr period before and after Coalinga is shown in (g). The radius of the Gaussian filter is 5 km, and the minimum number of shocks in the pre-Coalinga period needed to calculate a rate change, nmin, is 6. The seismicity rate cannot be reliably determined in areas rendered gray.

Within the first month after the Coalinga earthquake, the rate of microearthquakes increased in the creeping section and decreased along the Parkfield-Cholame section (Fig. 1b and Fig. 2d), except for a 2-week-long swarm that occurred 5 months after the Coalinga shock (at km=20 in Fig. 2d). During 1985-1990, the seismicity rate returned roughly to normal along the creeping section, but remained low along the Parkfield section (Fig. 2e), after which time the seismicity rate returned approximately to normal everywhere (Fig. 2f). To test whether these results are unduly influenced by the Gaussian smoothing, we also calculated the seismicity rate changes with only about 15% of the smoothing applied to Fig.2, as shown in Fig. 3. We reduced the radius of the filter from 5 km to 2 km; we also reduced n_min to 2, so that rates could be calculated in roughly the same areas as in Fig. 2. The same trends in the creeping (rate increases) and locked (rate decreases) section are evident, but there is more spatial variation.

Fig. 3 Calculated shear stress change (a), and observed seismicity rate change (b-d) on the San Andreas fault associated with the 1983 Coalinga-Nuñez shocks. Here the seismicity rate change has only ~15% of the smoothing used to produce Fig. 2, as the Gaussian radius is set to 2 km. The minimum number of shocks in the pre-Coalinga period needed to calculate a rate change, nmin, is lowered to 2. The seismicity rate change pattern is comparable to that shown in Fig. 2d-f, with rate decreases persisting with time at Parkfield, and rate increases diminishing with time in the creeping section. Although the functional dependence of seismicity rate change on stress change is higher than for the smoothed data (the slope, b, is 0.5), the correlation coefficient, R, is 0.32 between (a) and (b), and for the Coulomb stress change assuming m=0.4, it is 0.46, lower than for the smoothed data.

Several studies previously identified the seismicity rate decrease at Parkfield [Miller, 1996; Poley et al., 1987; Wyss et al., 1990], but the rate increase along the creeping section had escaped notice because it is briefer and spatially restricted. Wyss et al. [1990] described a seismicity rate decrease in the Parkfield region that began 2 years after the Coalinga shock. We find that the rate decrease southeast of Parkfield in fact began at the time of the Coalinga shock, but this was masked by the rate increase north of Parkfield, which did not end until about 18 months later. Thus the seismicity rate decrease that Wyss et al. ascribed to a process precursory to the next Parkfield earthquake we instead interpret as a postseismic response to the Coalinga shocks.

3. Calculated San Andreas Stress Changes

Several studies [Parsons et al., 1999; Reasenberg and Simpson, 1992; Stein, 1999; Toda et al., 1998] have found that seismicity rate change is correlated with the calculated Coulomb stress change, DCFF

DCFF = Dt + m Ds (1)

(where t is the shear stress and s is the normal stress, positive for unclamping). If correct, not only should aftershocks be more prevalent in regions of increased DCFF, but the rate of earthquakes should also drop in regions of decreased DCFF. The calculated shear stress change on the San Andreas fault imparted by the Coalinga-Nuñez events resembles the observed seismicity rate change along the fault at mid-crustal depths in map view (Fig. 1b) and in cross-section (compare Figs. 2b and 2d, or Fig. 3a and Fig. 3b). The correspondence suggests a causal relationship: San Andreas earthquakes became nearly twice as frequent (a log rate change of +0.3) where the stress increased by ~0.5 bars, and roughly half as frequent (log rate change of -0.3) where the stress decreased by the same amount. That the shear stress happened to have risen on the creeping section and dropped on the locked Parkfield-Cholame section (Fig. 1) is an accident of Coalinga’s location. Had the 1983 earthquakes struck north or south of Coalinga along any of the other frontal thrust faults, the inflection in stress applied to the San Andreas would not have coincided with the transition from the creeping to Parkfield sections.

The seismicity rate change and Coulomb stress change are statistically correlated, but the extent to which the correlation is driven by the shear stress change is not obvious. Visual inspection of Figs. 2b and 2d would suggest the shear stress controls the seismicity rate change (e.g., m~0 in eqn 1). However, a spatial regression of the Coulomb stress change on seismicity rate change (Fig. 4b) indicates that the regression coefficient, R, increases from m=0.0 to m=0.8; the regression for m=0.4 is shown in Fig. 4a. The functional dependence of the seismicity rate change on stress change (i.e., the slope, b), on the other hand, decreases from m=0.2 to m=0.8 (Fig. 4b). Thus the data lack the sensitivity to establish more than that 0.2²m²0.8. The same trends are evident in the data with minimal smoothing (Fig. 3), except that the slope b of the regression of seismicity rate change on stress change is higher, and the correlation coefficient R is lower (Fig. 4c). Regardless of the amount of smoothing or the assumed value of friction, the y-intercept is about -0.1, suggesting that the seismicity rates are biased toward negative values (i.e., in the absence of a stress change, the seismicity rates appear to decrease after May 1983). This negative bias is probably an artifact of the magnitude shift in the catalog that was not fully removed by using data starting from May 1980. Using only later does not circumvent this problem because the pre-Coalinga period would become too short to adequately measure seismicity rates.

Fig. 4 (a) Spatial regression of the calculated Coulomb stress change (with m=0.4), from Fig. 2c, on the observed seismicity rate change during the first 18 months after the Coalinga shock, from Fig. 2d. The number of earthquakes in 5-km-radius disks with centers spaced 1 km apart is counted and smoothed with a Gaussian function; thus each point is not independent of all others. Some 69% of the variance is accounted by an exponential fit to the data (black line, R=0.69). The positive Coulomb stress change data are also well fit by rate-state parameters appropriate for the creeping section, along which the stress changes are positive (gray curve). Dependence of the regression coefficients on friction, m, is shown for the highly smoothed (b) and lightly smoothed (c) data. The dependence of seismicity rate change on stress change (slope) is highest for small m, but the percentage of the variance explained by the regression (correlation, R) is highest for large m; thus these data do not permit an estimate of the friction coefficient.

4. Analysis of Surface Creep Changes

Surface creep slowed, stopped, or reversed on all creepmeters for 1-4 years after the Coalinga earthquake [Schulz et al., 1990], with the southern sites taking longest to recover (creepmeter locations are shown in Fig. 1a and Fig. 2a; the time series are shown in Fig. 5a). Fault creepmeter data permit independent analysis of the Coalinga stress transfer, and together with the seismic data, enable us to gauge parameters for a probability estimate.

Fig. 5 Response of creepmeters to the Coalinga earthquake, with records displayed north to south (top to bottom). (a) Observed and modeled surface creep; creepmeter locations are shown in Figs. 1a and 2. Rainfall is responsible for the rate increases in early 1983, but longterm creep rates and the creep retardation after Coalinga are judged reliable by Roeloffs (1999). A 0.6-mm offset occurred on the day of the 4 August 1985 Kettleman Hills shock on XMM1, and the creep rate reversed on XPK1 5-6 days later. No other changes are evident at the time of the Kettleman Hills shock. The end of the creep retardation is marked; the width of the marker reflects uncertainty in when the longterm creep rate resumes. The surface slip needed to shed the stress imposed by Coalinga, as calculated in Fig. 6b, is also indicated, with right-lateral slip positive. Subtracting the slip imposed by Coalinga from the longterm slip rate over the period of observed creep retardation produces the modeled creep series (dashed). (b) The linear relationship between the creep retardation and the longterm creep rate suggest a maximum retardation in the locked Parkfield section of ~4.4 yr; we equate this with the aftershock duration, ta, for the locked section.

The duration of retarded or reversed creep is correlated with the longterm creep rate (Fig. 5b), indicating that the faster the creep rate, the faster the recovery after the earthquake. A linear correlation, for which R=0.99, suggests that no retardation would occur where the fault creeps at the full slip rate of ~25 mm/yr, and that the recovery would last approximately 4-5 yr where the fault is fully locked, south of Cholame. A power law fit, in which the retardation period becomes infinite as the creep rate goes to zero, fits the data less well, with R = 0.84. An alternative is if the retardation period is instead a function not of the creep rate but the coseismic offsets. Although coseismic offsets of 0.1-1.8 mm accompanied the Coalinga shock [Mavko et al., 1985], the linear correlation between the retardation period and the observed offsets is poor (R = 0.05). Thus in what follows, we model the observed creep series by calculating the amount a frictionless San Andreas would slip in order to shed the stress imposed by the 1983 shocks, and subtract this induced slip from the longterm creep rate during the observed period of the creep retardation.

To find the distribution of San Andreas slip needed to relieve the stress imposed by the Coalinga earthquakes. We treat the crust as an elastic halfspace, and represent the San Andreas fault as a planar grid of boundary elements [Crouch and Starfield, 1983] that are free to slip except where the fault is locked. The locked region, where the fault slip rate is effectively zero, is approximately defined by an inversion of GPS data by Murray et al. [2000] (Fig. 6c). The fault patches used by Murray et al. are 2 x 2 km and extend to a maximum depth of 12 km. Such an elastic two-state model (where the fault is either free to slip or fully locked) nearly matches the decay of the longterm or secular creep rates toward the southeast (Fig. 6a). We then find the slip needed to shed the stress imposed solely by the 1983 Coalinga-Nuñez shocks (Fig. 6b). The effect of the 1983 shocks is to impose 6-25 mm of left-lateral slip, except at XSC1 where 9 mm of right-lateral slip is imposed. Finally, we subtract the imposed slip over the observed period of creep retardation (Fig. 5a). Thus creep reverses where the ratio of the imposed slip over the retardation period is larger than the longterm creep rate. At XMM1, the Coalinga shocks removed a calculated ~25 mm of San Andreas slip in one year; because 25 mm/yr is greater than the 16 mm/yr creep rate, the creepmeter moved left-laterally at 9 mm/yr for a year. The model matches most creep records; the largest departures are in the creeping section (XSC1), where accelerated creep is predicted but not observed, and at WRK1, where the observed retardation is larger than predicted.

Fig. 6 Calculated response of the San Andreas fault to secular loading (a) and to stresses imposed by the 1983 Coalinga-Nuñez earthquakes (b). The fault (c) is represented by a grid of freely slipping boundary elements except where it is fully locked, as inferred from inversion of GPS data by Murray et al. [2000]. The observed longterm creep rates (open circles fit by dashed line in a) are slightly lower than the model, which could be partly explained if the creepmeters do not span the full width of the fault zone. The southernmost sites are most sensitive to the detailed locking model. The calculated slip values in (b) are reproduced in Fig. 5a.

Finally, it is possible that some of the creep changes are the product of shaking associated with the Coalinga earthquake, rather than stress transfer. If shaking is responsible for the right-lateral coseismic jumps seen on all but XGH1, then the subsequent retardation could represent the time needed to return to the previously accumulated creep, with no long-term change in the creep trend. For all creepmeters except XSC1, however, there is a net decrease in the cumulative right-lateral creep since 1983 (Fig. 5b). These creep deficits are close to the values calculated in the boundary element model (Fig. 6b), from 25 mm at XMM1 to 6 mm at XHG1. Thus while it is difficult to eliminate the possibility of contamination by shaking in these records, the data are more consistent with stress interaction.

5. Comparison with Other Studies of the Creep Series

Simpson et al. [1988] fit the creepmeter records to a model that is, like ours, driven by stressing from the deeper San Andreas fault, and modulated by the stress change associated with the Coalinga earthquake. In their model, the San Andreas stressing rate below 5-km depth is 0.75 bar/yr at XMM1 and 0.5 bar/yr at XGH1. They let the upper 5 km of the fault respond in a linear viscous manner to the stress changes, producing a left-lateral excursion at XMM1, and smaller left-lateral excursions on creepmeters to the south, all of about one-year duration. There is no tendency for longer creep retardation periods toward the south, but the magnitudes of the creep changes resemble the observations. It is unclear, however, whether their modeled creep reversal occurred because they doubled the Coalinga coseismic slip to make the changes in the creep rates more apparent.

Examining the 1983 coseismic offsets in the creep series, Mavko et al. [1985] found a rough match between the creepmeter offsets and calculated Coulomb stress changes for the Coalinga-Nuñez shocks, using m=0.6. Recent analysis by Roeloffs [2001] of the response of creepmeters to rainfall and earthquakes indicates, however, that re-centering of the instruments during shaking contaminates the coseismic displacements, whereas the creep rates before or after earthquakes suffer fewer such problems.

6. Calculation of Earthquake Probability

The time-dependent response of seismicity and creep to the stress imposed by the Coalinga-Nuñez sequence can be incorporated into an earthquake probability calculation [Stein, 1999]. This approach differs from a statistical analysis of the probability that a small (2<M<6) Parkfield shock is likely to followed by a M~6 mainshock [Michael and Jones, 1998]. In Coulomb failure theory, a positive or negative stress change on the San Andreas fault causes an advance or delay to the time until failure is reached (Fig. 7, upper panels). This results in a modest but permanent change in earthquake probability (Fig. 7, thin solid lines in lower panels). If this were a complete description of the process, then the Parkfield earthquake would be delayed by the stress change (-0.15 bars) divided by the stressing rate (~0.1 bar/yr), or 1-2 yr. While reasonable, such an approach fails to explain seismic observations for several other faults. Examples include the order-of-magnitude decrease in M>6 seismicity in the San Francisco Bay area during the 75 years after the great 1906 earthquake [Bakun, 1999; Harris and Simpson, 1998], and the 12 progressive M>6.7 earthquakes along 1000 km of the North Anatolian fault in 60 years starting in 1939 [Barka, 1996; Stein et al., 1997]. In such cases, the calculated static stress changes of several bars would at most advance or delay subsequent earthquakes by decades, and could not explain such persistent seismicity decreases or increases persisting for 60-75 years.

Fig. 7 Schematic illustration of the effect of negative (at point x) and positive (at point y) stress changes on earthquake probability along a strike-slip fault, given an arbitrary off-fault earthquake source (rupture). The 10-yr probability is the integral over the period in question. Rate and state friction effects associated with the shift to a time earlier or later in the earthquake cycle are ignored in this example. The lower panels can be compared to the observed seismicity rate changes in both the creeping and Parkfield sections (Fig. 9) and calculated probabilities for Parkfield (Fig. 10).

One solution is incorporation of rate and state friction into the probability model. In rate/state friction, seismicity is viewed as a sequence of nucleation events in which the state depends on the fault slip, slip rate, and elapsed time since the last event [Dieterich, 1994; Dieterich and Kilgore, 1996]. The seismicity rate equation is

(2)

The transient effect of a stress decrease strongly amplifies the permanent decrease [Toda et al., 1998], because the fault slips at a lower rate, causing a lower rate of earthquake nucleation (Fig. 7, lower panels). The transient recovery time is inversely proportional to the fault stressing rate times As, where A is a constitutive parameter, and s is the total normal stress [Dieterich, 1994; Dieterich and Kilgore, 1996]. In addition to the stress change, to evaluate (2) one must estimate the aftershock duration or recovery time, and the fault stressing rate or As. For an interaction-based renewal probability, one further assumes that with the passage of time from the last M~6 shock in 1966, another such earthquake becomes more likely. To calculate an interaction-based renewal probability, one must assume a probability density function and also estimate the elapsed time since the last earthquake, an inter-event time, and the coefficient of variation for such events. Despite inevitable uncertainty in such assignments, the short repeat time, similar size, and long historical record for Parkfield earthquakes make estimates more reliable here than for most faults.

Transient decay. The decay is proportional to the aftershock duration, t_a, the time elapsed until the rate of seismicity returns to the rate that prevailed before the mainshock occurred. In rate and state friction, the aftershock duration t_a is independent of mainshock magnitude, and t_a is instead related to As through

t_a =As /(3)

where is the fault stressing rate [Dieterich and Kilgore, 1996]. We use the four largest mainshocks on the San Andreas in the NCSN catalog to estimate aftershock duration. We calculate the aftershock decay rate for M>1.3 shocks, and use the M>1.3 seismicity rate from 1980 until the time of mainshock in the same zone to estimate the background rate (Fig. 8a-c). For the 1966 mainshock (Fig. 8d), we use the catalog of Meagher & Weaver [2000], and measure the background rate starting in 1937 for M>3.5, the completeness level of this catalog. There is an apparent increase in aftershock duration to the southeast, with aftershock durations growing from ~0.6 yr in the creeping section to ~5 yr in the locked region. Because we are limited by the smaller maximum size of earthquakes in the creeping section, it is also possible that the aftershock duration is a function of magnitude. We note, however, that the measured aftershock durations are comparable to the observed creep retardation periods (shown in the map panel in Fig. 8), suggesting that the seismic and creep observations are manifestations of the same process of transient recovery to sudden stress changes. Thus as mean values, we will take t_a to be ~0.5 yr in the creeping section and ~4.0 yr in the locked Parkfield section.

Fig. 8 Estimated aftershock durations, ta, along the San Andreas fault, with plot features labeled in (a) and earthquake and creepmeter locations shown in the map panel at bottom. The largest earthquakes along the creeping and Parkfield sections were selected at the four sites to permit the best possible calculation of the aftershock decay. Aftershocks were extracted from the NCSN catalog in the boxed regions to produce the rate decays; the background seismicity rate was estimated from the shocks in the same boxes from 1980 until the time of the mainshock. For the 1966 Parkfield earthquake, the catalog of Meagher and Weaver [2000] was used, with the background rate estimated for 1937.5-66.5; and from the NCSN catalog for 1980-2000. There is a general trend toward longer durations in the more fully locked part of the fault, consistent with the creep retardation periods from Fig. 5a.

Although far from proven, these t_a assignments are consistent with another means to gauge the aftershock duration. If the stressing rate, , can be approximated by the mainshock stress drop, Dt, divided by the inter-event time, t_r, then from (3)

t_a = t_r (As/Dt) (4)

In the Parkfield section, the maximum observed earthquake magnitude M_max~6 and t_r~22 yr. Nadeau and Johnson [1998] find that in the creeping section, M_max~4 and t_r~2.4 yr. Thus if earthquake stress drops and As are constant in both locations, one would expect t_a to be roughly an order of magnitude larger in the locked section, consistent both with the observed creep retardation periods and aftershock durations.

Fault stressing rate. In order to calculate the interaction-based probability, we need to estimate the fault shear-stressing rate , shown schematically in the middle panels of Fig. 7, at sites where earthquakes occur. Three conceptual approaches to estimating lead to different answers. If one assumes that the San Andreas is vertical, straight and subject to the same plate boundary tractions–or alternatively, the same deep slip rate and locking depth–throughout the region of Fig. 1, then its tectonic shear stressing rate should be uniform along strike. At a mid-crustal depth of about 8 km, the shear stressing rate would be about 0.1 bar/yr. If one instead assumes that where the fault creeps, stress is relieved, then the stressing rate would be near-zero north of XMM1 and much higher than 0.1 bar/yr at the north end of the locked zone, where dislocations would continuously pile up (Fig. 6c). Another alternative, inspired by the work of Rubin et al. [1999] and Waldhauser et al. [1999], is that seismicity in the creeping section is concentrated along isolated streaks, with the patches between streaks undergoing steady, aseismic creep. If this were correct, then the stressing rate in the streaks might be higher than in the Parkfield locked patch, where the load could be more uniformly distributed. (Earthquake streaks are not evident in Fig. 2d-g because pre-1984 data cannot be relocated by waveform cross-correlation.)

Because we can not confidently eliminate two of these alternatives, we instead find the stressing rate that matches the observed seismicity rate as a function of time, and satisfies the creep retardation and aftershock duration data. The creeping section (Fig. 9a) sustained a calculated mean 0.3-bar Coulomb stress increase, and t_a there is 0.5-1.0 yr on the basis of Fig. 8. The decay of the seismicity rate as a function of time is best fit in the state/rate formulation with a stressing rate of 0.5 bar/yr, about five times higher than the average San Andreas stressing rate, and thus a rate appropriate for seismic streaks. This stressing rate is also compatible with the spatial regression of seismicity rate on stress change shown in Fig. 4 (gray curve). Thus both the temporal and spatial seismicity rate data are well fit with a high stressing rate in the creeping section. In the Parkfield section (Fig. 9b), the mean calculated Coulomb stress change is —0.15 bars, and the data are satisfied by a stressing rate of 0.1 bar/yr for t_a = 2-4 yr. Thus, given one degree of freedom (the stressing rate), we can satisfy the temporal decay of the seismicity rate following the Coalinga-Nuñez earthquakes.

Fig. 9 Observed and theoretical seismicity rate as a function of time, showing the response of the San Andreas fault to the 1983 Coalinga-Nuñez earthquakes. The observed seismicity rate changes are averages for the rectangles shown in Fig. 2d-f and Fig 3b-d.These values have been shifted by a log rate change of —0.1, because of the negative seismicity rate bias discussed in the text, and indicated in Fig. 4 as the Y intercept. (a) The mean Coulomb stress change (DCFF) in the creeping section, from Fig. 2c, is 0.3 bar. The modeled transient change in seismicity rate as a function of aftershock duration, ta (for ta=0.1 yr, As=0.05 bar; for ta=0.5 yr, As=0.25 bar; for ta=1.0 yr, As=0.50 bar. (b) The mean Coulomb stress change in the Parkfield section is —0.15 bars (for ta=1 yr, As=0.1 bar; for ta=2 yr, As=0.2 bar; for ta=4 yr, As=0.4 bar). The first interval covers 0.5-1.0 yr after the Coalinga shock to avoid the time period of the swarm at km=22 in Fig. 2d. The seismicity rate change data are approximately fit by ta~2-4 yr, in accord with creep (Fig. 5b), and seismic data for the 1966 aftershock duration (Fig. 8d).

Parkfield probability. The probability follows naturally from the seismicity rate change plot of Fig. 9. We assume that the seismicity rate changes not just for microearthquakes but at all magnitudes. Because M~6 earthquakes are infrequent, there is only a chance that the rate change will result in a detectable change in the occurrence of a M~6 event after 1983. We thus perform a Monte Carlo analysis of 1,000 runs in which the tested values are drawn from a Gaussian distribution of the input parameters, and plot the mean value and uncertainty as a function of time (Fig. 10). The calculated stress change at the 6 x 6 km hypocentral site of the 1934 and 1966 earthquakes is —0.3±0.1 bar (Fig. 2b-c). The fault stressing rate (0.1±0.025 bar/yr), aftershock duration (4±1 yr), and As (0.4±0.1 bar) are estimated from the preceding analysis of the creep and seismicity data. Calculations are made alternating between lognormal and Brownian Passage Time [Matthews et al., 2001] probability density functions, with mean inter-event time of 22 years, and a coefficient of variation of 0.5 (given a range of 0.35 [Savage, 1993] to 0.70 [Roeloffs and Langbein, 1994]). The calculated 10-year probability (Fig. 10) in 1983 decreased from 54±22% to 42±23%, and is not calculated to have returned to the pre-Coalinga probability until about 1991. October 1992 marked the beginning of a period of heightened seismic activity at Parkfield [Michael and Jones, 1998] (including a M=4.7 shock in 1992, M=4.4 and M=4.8 shocks in 1993, and a M=5.0 shock in 1994), the occurrence of which is consistent with our probability calculation. Although the probability undergoes a fractional drop of 22% in 1983, it is not statistically significant. The rate/state effects of the stress change could be added to a Poisson model, as was done by Toda et al. [1998]. This would yield smaller nominal uncertainties because the coefficient of variation of the inter-event time (0.5) would not enter into the calculations. However, we regard renewal as a better description of earthquake recurrence at Parkfield.

Fig. 10 Ten-year probability of M>6 earthquakes at Parkfield, calculated by Monte Carlo analysis following Parsons et al. [2000]. The 22-yr mean inter-event time for M~6 shocks is what one would have adopted in 1983 (since the earthquake has not struck for 34 yr, one would use a slightly longer inter-event time today). The 10-yr probability drops from 54±22% to 42±23% as a result of the stress decrease imposed by the Coalinga earthquake, and does not recover to pre-1983 values of probability until 1991. The 10-yr probability of such an earthquake for the ensuing decade, 2001-2011, is calculated to be 58±17%.

7. Conclusions

The thrust of our findings is that stress increases and decreases associated with nearby earthquakes influence subsequent seismicity. Seismicity rates are easily observed at Parkfield; they change at the time of the 1983 Coalinga-Nuñez shocks in a manner that resembles the calculated Coulomb stress change imparted by the 1983 shocks. Although we can explain the observations without a viscous[Ben-Zion et al., 1993; Freed and Lin, 1998; Simpson et al., 1988] or poroelastic rheology [Miller, 1996], such processes may indeed be important. Nevertheless, a probability model governed by rate and state friction and driven by steady stress buildup and stress transfer from nearby earthquakes can satisfy most of the Coalinga-Parkfield observations. We find the effect of the stress decrease on earthquake probability at Parkfield was larger and lasted longer than previously supposed. Whether the absence of a Parkfield earthquake since 1983 is the result of the Coalinga-Nuñez events will never be known, but we calculate that the probability of a M~6 Parkfield shock today is higher than it was before 1983.

Acknowledgements. We are indebted to Robert Simpson, Paul Reasenberg, Felix Waldhauser, Stefan Weimer, Tom Parsons and Evelyn Roeloffs for their analytical tools, insight, and reviews. We gratefully acknowledge funding from PG&E.

References