Space–Time Asymmetry in Earthquake Pairs along the Central San Andreas Fault: Evidence for Small Earthquake Links at Long Distances

California’s central San Andreas fault (CSAF), between Parkfield to the southeast and San Juan Bautista to the northwest, is a well‐studied fault segment that is distinguished by ongoing and substantial shallow fault creep. The creep rate varies along strike, reaching the maximum of about 25–35 mm/yr (e.g., Ryder and Bürgmann, 2008; Maurer and Johnson, 2014; Scott et al., 2020; Liu et al., 2022) south of Bitterwater, where it nearly equals the long‐term fault slip rate and decreasing to the north toward San Juan Bautista. Although the long‐term average creep rate has been nearly constant since 1970 (Ryder and Bürgmann, 2008), creepmeters have detected numerous individual creep events with durations of hours to days (e.g., Gittins and Hawthorne, 2022). InSAR analyses have also suggested temporal changes in creep rate at time scales of months to years (Khoshmanesh et al., 2015; Khoshmanesh and Shirzaei, 2018), some of which appear to have a seasonal component (Li et al., 2022). The CSAF is characterized by active seismicity with over 40,000 earthquakes recorded by the Northern California Seismic Network (NCSN) since 1981. Most of these events occur between 2 and 8 km depth. Thurber and Sessions (1998) found a statistical correlation between creepmeter recorded events and M > 3.8 earthquakes in the ensuing 5–10 days. Nadeau and McEvilly (2004) identified clusters of repeating earthquakes along the CSAF with variations in repeat times that suggest quasi‐periodic variations in deep slip rate, assuming a model in which repeating earthquakes are driven by adjacent aseismic regions with slip that repeatedly loads the rupture region to failure. Li et al. (2022) found a correlation between changes in deep slip rate along the CSAF inferred from similar earthquake repeat times and InSAR measurements of surface creep rate.

Here, we search for pairs of earthquakes along the CSAF over four decades that occur within a minute of each other. Most of these nearly contemporaneous events are physically separated by less than a few kilometers. However, when we examine the small number of events with spatial separations of more than 5 km, we identify a surprising asymmetry. For event pairs with temporal separations of 5–50 s and spatial separations of 5–50 km, it is much more likely that the later event will be northwest of the earlier event than it will occur to the southwest. This asymmetry is not predicted by standard earthquake‐to‐earthquake triggering models and appears unrelated to any dynamic triggering from distant events. We explore possible explanations for these observations but are unable to identify a definitive cause.

We use earthquake locations and origin times from the Northern California Earthquake Catalog as downloaded from the Northern California Earthquake Data Center (NCEDC, 2014). Figure 1 shows the location of our study region along the CSAF, which contains 45,868 earthquakes between 1981 and 2021. The along‐strike distance along the SAF is defined from a southeast reference point at (36.20° N, −120.75° E) toward a northwest reference point at (36.90° N, −121.64° E). We include all earthquakes within ±5 km from this reference line.

Temporal clustering of seismicity is commonly observed, and has origins in both aftershock triggering and swarm activity. Here, we are interested in pairs of earthquakes occurring close in time that are not part of the immediate foreshock or aftershock sequence of a large earthquake (e.g., Meng and Fan, 2021). Accordingly, we first remove from the catalog any earthquakes within the 0.5 day before or the 3 days after the five M ≥ 6 earthquakes in the regional catalog and the two M ≥ 5 earthquakes within our CSAF study region. Next, we search for earthquake pairs in the CSAF that are within 150 s in time and within 100 km in epicentral distance, such that only horizontal distance is used and any depth difference is ignored.

Figure 2a plots the time and distance offset of the second event in each pair with respect to the first. Positive distances are those to the northwest, and negative distances are those to the southeast. The vast majority of these earthquake pairs are separated by less than 5 km. However, we are interested in the smaller number of earthquake pairs that occur at larger spatial offsets. These event pairs exhibit a pronounced time–distance asymmetry in the first minute, with more pairs in which the later event occurs to the northwest than those in which the second event occurs to the southeast. The asymmetry is particularly pronounced within a 5–50 s and 5–50 km window, shown as the rectangles in Figure 2a. Within this window, 42 of the 51 event pairs involve northwestern directivity.

Figure 2b plots catalog magnitude versus time for these 51 earthquake pairs, compared to the other events in the NCSN catalog. Blue is used for the northwest moving event pairs and red for the southeast moving event pairs. The earthquake pairs occur throughout the time period, with most magnitudes between about 0.5 and 2.5. Interestingly, the first southeast offset pair does not occur until 1997, following 23 consecutive northwest offset pairs starting in 1981. Two pairs of events on 3 May 1995 have the same second event, with their initial events occurring about 20 s and less than one kilometer apart.

Next, we explore the spatial distribution of these earthquake pairs compared to the CSAF seismicity. Figure 3 shows our CSAF study region in map view and depth cross section, with the event pairs connected with colored lines in the depth cross section. Northwestward offset event pairs are shown in blue and southeastward offset pairs in red (Fig. 3b,c). Observe the dominance of the northwest offset pairs (blue), especially in the years before 1997 and northwestward distances less than 70 km.

Taken at face value, the along‐strike asymmetry observed for the earthquake pairs within our time–distance window appears unlikely to have occurred due to random change. Assuming that each direction has equal probability, the chances that 42 out of 51 pairs would be in the same direction (either northwest or southeast) is about 0.00034% (about 1 chance in 300,000). This result is robust with respect to our aftershock declustering that removes seismicity within a few days of the local M > 5 earthquakes. If we include all earthquakes in the catalog (i.e., do not remove any aftershocks), we find that 43 out of 54 pairs have northwest offsets, which would occur (in either direction) with only 0.0014% probability by random chance (about 1 chance in 71,000). If we remove more possible aftershocks by reducing our mainshock minimum magnitude from 5 to 4.5 within our study region, we find that 38 out of 47 pairs have northwest offsets, which would occur (in either direction) with only 0.0025% probability by random chance (about 1 chance in 40,000).

However, these calculations ignore the “data fishing” aspect of our choice of this fault segment and the space–time window boundaries in Figure 2. We have searched for near‐contemporaneous but well‐separated (>5 km) event‐pair asymmetry for five or so other fault segments in California with substantial ongoing seismicity and obtained suggestive but inconclusive results in some cases. Because the CSAF stood out as being particularly anomalous, we have focused our attention on trying to understand its behavior first.

How important are the exact locations of the box boundaries shown in Figure 2a for our statistics? This can be evaluated by craftily adjusting the limits to just barely include or not include events. As shown in Figure S1, available in the supplemental material to this article, one can make the observations appear more statistically significant by changing the box boundaries slightly to 5.25–50 km and 5–53 s (shown as the dashed blue box), which will enclose 46 blue pairs out of 54 total, which would occur randomly (in either direction) just 1 out of 7.2 million times. Alternatively, somewhat larger changes of the box boundaries to (5–38 km, 13–63 s, shown as the dashed red box) will enclose 35 blue pairs out of 51 total, which would occur randomly 1.09% of the time. This is a huge difference, but the last result remains statistically significant using standard thresholds (e.g., the 95th percentile). What about larger adjustments in our windowing choice? The asymmetry disappears at times greater than about 60 s, so including some of these points would reduce the statistical significance of our results. As we will discuss later, this lack of asymmetry at later times argues against slow‐slip explanations for our results, because slow‐slip events are observed to propagate at much slower velocities than those associated with our space–time window.

At closer distances, the asymmetry is still present, in that 58.5% of the 53 pairs within (2–5 km, 5–50 s) involve northward propagation. However, the lesser strength of the asymmetry means that this would have occurred by random chance about 14% of the time, that is, it would not be considered statistically significant by itself. Including these closer events together with more distant events reduces the significance of the asymmetry compared to the more distant window alone, but the result remains very unlikely to be due to random chance. A window of (2–50 km, 5–50 s) encloses 73 blue pairs out of 104 total, which would occur randomly (in either direction) only 0.047% of the time (1 chance in 21,000). However, it has always been the case that we have been searching for correlations among pairs separated by at least 5 km, because we are interested in possible distant earthquake interactions that are not explained by near‐source aftershocks or location error. We did not exclude the pairs separated by less than 5 km at a later stage in the analysis to increase the statistical significance of our results.

From a Bayesian perspective, probability estimates should also take into account the plausibility of a result and the lack of a physical model to explain our observed asymmetry suggests that a high standard of proof should be applied. Taken all of these factors into account, it is difficult to objectively assign a precise level of statistical significance to our results for the CSAF, and we leave it to the reader to judge how likely the observed asymmetry is caused by random chance. The ultimate significance test will be to see what happens in the future within our time–distance window (5–50 km, 5–50 s), but it likely will take many years to obtain a definitive result, given the slow rate of earthquake pair occurrence on the CSAF.

Assuming the observed asymmetry is real and not due to random chance, what could be causing this signal? One possibility is that it reflects a bias in the earthquake catalog, possibly caused by nonuniform station coverage. For example, southward‐propagating event pairs may be easier to detect at stations to the north, because their arrivals are more spread out in time (so the second‐event arrival is further from the coda of the first event) and harder to detect at stations to the south. Thus, if there were more stations to the south, this might provide an explanation for the asymmetry. To investigate this possibility, we downloaded phase‐pick information (NCEDC, 2014), which gives the stations that were used to locate each of the events in our earthquake pairs. Plots of the earthquake and station locations for the 51 pairs are included in Figures S3–S7.

These plots show that the number of contributing stations increases with event magnitude, as might be expected. However, overall there are somewhat more contributing stations to the northwest than to the southeast, which would be more rather than less likely to detect southeast‐propagating pairs. The smaller events are generally located using stations within about 30 km around them, and there is no obvious network asymmetry that could cause a preference for northwest‐propagating pairs. For example, the northwest‐propagating 22 April 1990 pair (Fig. S3, bottom right) appears equally likely to have been recorded if the event order was switched to make it southeast‐propagating, and the southeast‐propagating pair on 31 March 2017 has similar station coverage to the northwest‐propagating pair on 5 April 2018 (see Fig. S6). We also plotted the number of NCSN picks versus magnitude (Fig. S2) and see no significant difference between the first‐ and second‐event pick numbers, suggesting that the coda of the first event does not noticeably limit the number of picks for the second event, and we also see no clear difference in second‐event pick numbers between the northwest‐ and southeast‐propagating pairs. Overall, we find no evidence that network geometry is causing the observed pair asymmetry, but further analysis of possible catalog biases may be warranted, given the improbability of other explanations for the asymmetry (see the following).

Earthquake occurrence is more clustered in time and space than purely random catalogs predict. This clustering results from both earthquake‐to‐earthquake triggering (i.e., aftershock sequences and ETAS‐like models as well as long‐range dynamic triggering) and physical driving mechanisms such as the fluid flow or slow slip that are thought to drive many earthquake swarms. However, most seismicity clustering models do not predict long‐range correlations in earthquake occurrence or the along‐fault pair asymmetry that we observe (i.e., that the later events in the pairs are more likely to be northwest of the initial event).

The magnitudes of our earthquake pairs are small enough (see Fig. 2) that one would not expect the first event to generate many aftershocks, particularly at distances of 5–50 km. If earthquake‐to‐earthquake triggering is involved, it seems more likely that this is caused by distant earthquakes, in which case the asymmetry might result from earthquakes to the southeast that generate northwest‐propagating body or surface waves that sequentially trigger the events. However, Hirao et al. (2021) found no clear earthquake triggering in the CSAF area resulting from the 2003 $Mw$ 6.5 San Simeon and 2004 $Mw$ 6.0 Parkfield earthquakes to the south. The short temporal separations between the first and second events could have been due to passing seismic waves from remote earthquakes, which may have near‐simultaneously dynamically triggered the two events. To explore this hypothesis, we examine regional and teleseismic catalogs at the times of our 51 event pairs. For each event pair, we search the International Seismological Centre catalog (International Seismological Centre, 2022) to identify the closest M ≥ 6 global earthquake preceding the pair. The epicentral distance to the mean location of the event pair is then examined with the preceding time (Fig. 4). Assuming an average surface‐wave velocity of 3.5 km/s, we find no correlation with predicted ground motions from large earthquakes. None of the events in the pairs occurs within one hour following the surface‐wave arrival of an M ≥ 6 earthquake occurring anywhere on the globe (see Fig. 4). In addition, the time versus distance offset of the event pairs is generally much slower than surface‐wave velocities (see Fig. 2a). The findings suggest that the event pairs are unlikely instantaneous triggering caused by transient external perturbations due to distant earthquakes.

However, the seismic waves may have triggered creep that triggered earthquakes at later times (e.g., Shelly et al., 2011). Three of our observed pairs occurred within one day of the passing waves of a distant M ≥ 6 global earthquake, and 28 pairs occurred within one week afterward (Fig. 4). Thus, triggered creep is a possible scenario that could explain the delay between the seismic‐wave passage and the occurrence of the earthquake pairs. However, as shown in Figure 4 (inset), the interevent times of sequential global M ≥ 6 earthquakes center around days to a week, suggesting that the fact that 28 pairs occurred within one week of an M 6 earthquake is more likely due to random chance.

We find the significant seismicity rate episodes tend to cluster in time, and the four segments are not always synchronized (Fig. 5b–e). Interestingly, the significant seismicity rate episodes do not always coincide with local M ≥ 4 earthquakes with 74% of the earthquakes accompanied with three‐day seismicity anomalies. In comparison, 25% of the 51 event pairs are associated with three‐day seismicity anomalies. The correlation suggests that some of our observed earthquake pairs may have been part of clusters that involve additional earthquakes. However, these clusters do not seem to have clear spatiotemporal migration patterns. Moreover, these clusters lack obvious mainshocks and do not exhibit Omori‐law decay from the initiating earthquakes. The clusters may have been swarms, which are typically thought to be driven by fluid flow or slow slip and often exhibit spatial migration. However, both of these mechanisms typically do not operate fast enough to trigger a second event that is tens of kilometers away from the first event within 50 s.

Geodetic data suggest that slow‐slip events typically have along‐strike propagation velocities of 5–15 km/day (Behr and Bürgmann, 2021). At Parkfield (along the SAF south of Bitterwater), strainmeter and creepmeter observations suggest slow‐slip propagation velocities of 4–45 km/hr (Mencin, 2018). Along‐strike migration has been observed in tremor, and families of deep low‐frequency earthquakes (LFEs) observed along the SAF south of Bitterwater (e.g., Shelly, 2015, 2017). Some of this LFE activity has been correlated to slow‐slip events observed with borehole strainmeters (Delbridge et al., 2020). However, these LFE events are south of and occur at much greater depths (15–30 km) than our earthquake pair observations (mostly at 5–10 km), so their relevance is unclear. In addition, observed LFE propagation velocities are typically 25–100 km/hr, much slower than the apparent velocities of 0.2–4 km/s consistent with our earthquake pairs.

To check if there are any observed strain anomalies associated with our earthquake pairs, we searched borehole strainmeter records at stations B058, B065, B066, and B067 (see Fig. 3 for locations) between 2008 and 2019, to see if there are any obvious signals at the times of the event pairs that might suggest the occurrence of slow‐slip events. However, we find no identifiable signals, suggesting that any slow‐slip events at the time of the pairs are too small to observe with the strainmeters. We have also examined the CSAF creep‐event catalog of Gittins and Hawthorne (2022), and find no clear correlation between observed creep events and the times of our earthquake pairs. In addition, observed creep‐event propagation velocities are all less than 500 km/day, with most velocities less than 100 km/day.

It is interesting that observations of microearthquake aftershocks located within a few hundred meters of their mainshocks along the CSAF also suggest preferential occurrence to the northwest (Rubin and Gillard, 2000; Rubin, 2002; Rubin and Ampuero, 2007). This aftershock asymmetry has been hypothesized to result from stress asymmetry caused by preferential rupture propagation to the southeast, which may result from a bimaterial interface with a velocity contrast across the fault. However, these stress changes would not extend far enough to explain the more distant (5–50 km) earthquake pair asymmetry shown in Figure 3.

The microearthquake aftershock studies used waveform cross correlation to provide precise relative event locations among nearby events (Rubin and Gillard, 2000; Rubin, 2002; Rubin and Ampuero, 2007). Such relocated catalogs are necessary to provide enough location accuracy to resolve details at the scale of tens to hundreds of meters. There is a double‐difference earthquake catalog for northern California (Waldhauser and Schaff, 2008; Waldhauser, 2009) that could be used for this purpose. However, it currently only goes back to 1984, and thus we prefer to use the standard NCSN catalog to provide uniform coverage during the entire 1981–2021 time period. Because we are examining earthquake pairs with at least 5 km separation, precise location accuracy is not required.

Given the lack of an obvious explanation, we now indulge in some speculation regarding what could cause the earthquake pair space–time asymmetry, recognizing that these hypotheses lack any direct supporting evidence and may seem physically implausible. Two ideas are:

1. The pairs are triggered by a new class of aseismic slip events that propagate at near‐seismic velocities. In this case, the observed asymmetry might result if these events tend to originate in the more freely slipping SAF segment to the south of Bitterwater (where the creep rate nearly equals the long‐term slip rate) and then propagate into the partially coupled segment of the SAF north of Bitterwater. Alternatively, these slip events might have a preferred propagation direction due to across‐fault velocity contrasts, although it should be observed that earthquake rupture directivity related to across‐fault velocity contrasts (e.g., Andrews and Ben‐Zion, 1997) is caused by dynamic effects for slip pulses, which presumably do not apply to aseismic events.

2. The second event in many of the pairs is dynamically triggered by the first event. Dynamic triggering could explain the small time separation between the events, which is suggestive of a fast apparent propagation velocity, but the generally small magnitudes of the first event and large distances to the second event pose challenges to this idea. Conceivably dynamic triggering might be encouraged by fault‐zone‐guided waves excited by the first event, which would experience less amplitude decay with distance than crustal P and S waves. In this case, the observed asymmetry might result if CSAF earthquakes have a preferred northwest rupture direction (i.e., opposite to that proposed by Rubin and Ampuero [2007] to explain aftershock patterns) leading to higher radiated amplitudes to the northwest because of directivity. Alternatively, the CSAF low‐velocity‐zone structure might vary along strike in such a way that it focuses or preserves energy more efficiently for northward wave propagation than for southward propagation.

In the future work, we plan to explore whether the earthquake pair asymmetry that we observe is unique to the CSAF, which is unusual for its ongoing fault creep, or whether other faults also exhibit this behavior. It should be observed that the earthquake pairs that best show the asymmetric signal are of small‐magnitude events and only occur about once per year. Thus our study benefited from the large number (over 45,000) of earthquakes detected and located by the NCSN within our study region over more than four decades.

The earthquake catalog for this study was accessed through the Northern California Earthquake Data Center (NCEDC; doi: 10.7932/NCEDC). The supplemental material to this article includes Figures S1–S7.

The authors thank Roland Bürgmann and David Shelly for their insightful and constructive reviews. Funding for this research was provided by U.S. Geological Survey (USGS) Grant G20AP0027 and by the Southern California Earthquake Center (SCEC) under Grants 20097 and 21055. SCEC is funded by NSF Cooperative Agreement EAR‐1033462 and USGS Cooperative Agreement G12AC20038. Wenyuan Fan was supported by NSF‐EAR 2022441. This is SCEC contribution number 12693.