of 24
Supplementary Materials for
A unified perspective of seismicity and fault coupling along the
San Andreas Fault
Yuan-Kai Liu*, Zachary E. Ross, Elizabeth S. Cochran, Nadia Lapusta
*Corresponding author. Email: ykliu@caltech.edu
Published 23 February 2022,
Sci. Adv.
8
, eabk1167 (2022)
DOI: 10.1126/sciadv.abk1167
This PDF file includes:
Materials and Methods
Supplementary Text
Figs. S1 to S8
References
Materials and Methods
Data:
The dataset used in this work is the Northern California Seismic Network double-difference
(DD) earthquake catalog (
15, 56, 57
). The automatic and routinely relocated earthquake cata-
log was downloaded from the Northern California Earthquake Data Center. The catalog covers
events of the entire central and northern California since January of 1984 to February of 2021.
We further select a subset of this catalog along the strike of the 225-km-long central San An-
dreas Fault (central SAF) extending from around Loma Prieta in the northwest (122.25
W,
37.40N) to near the Carrizo segment in the southeast (120.14
W, 35.59
N) within
±
5
km dis-
tance from the fault. The majority of the hypocenters locate at around 2 to 8 km depth (Fig.
S1). To include the deeper events from the 1989 Loma Prieta sequence which sit on a dipping
fault plane, we increase the width of our data selection to
±
10
km for the region north of San
Juan Bautista.
The regions further north of the Loma Prieta and further south of Cholame (the Carrizo
segment) are dominated by rarely occurring historic large mainshocks and aftershocks. These
include the 1857 M7.9 Fort Tejon earthquake (
58
) in the south, 1906 M7.9 San Francisco earth-
quake (
59
) in the north. The above segments at the both ends of our study profile have insuf-
ficient events during periods devoid of large shocks for interevent-time analysis and thus can
easily lead to overestimation of the fraction of non-clustered seismicity due to the sampling bias
(
32
). We hereby exclude these portions of the central SAF from our discussion.
Magnitude of completeness estimation:
The magnitude of the catalog is duration magnitude (
M
d
). We assess the magnitude of
completeness (
M
c
) using the maximum curvature (MAXC) method (
60
). This method defines
the point of the maximum curvature as the magnitude of completeness by computing the max-
imum value of the first derivative of the frequency-magnitude distribution curve. In practice, it
is equivalent to the magnitude bin with the highest frequency of events in the non-cumulative
frequency-magnitude distribution. Since the MAXC technique tends to slightly underestimate
the magnitude of completeness, we add 0.2 for correction (
61
). Based on the selected subset
of the catalog, we estimate
M
c
for every year from 1984 to 2021 and find that the values are
within the range of 0.8–1.5 (Fig. S2). We thus set a constant minimum magnitude
M
min
=1.5
for the entire analysis, which results in 24,082 events over 37 years (Fig. S3). With this subset
catalog, we apply a moving spatial window along the fault (15 km along-strike by 10 km fault-
normal, step size 1 km) and compute properties of seismicity within the spatial window. We
also make sure that, within each window, there are more than 1,000 events larger than
M
c
. The
computed seismicity properties, including b-values, non-clustered seismicity rate, and fraction
of non-clustered seismicity, are described in the following sections.
b-value estimation:
We estimate the b-value of the dataset with the maximum likelihood method proposed by
Aki (
62
), where, for each subset of events within the spatial window, the b-value is computed
as:
b
=
log
10
(
e
)
M
M
min
where e is the Euler’s number,
M
is the mean magnitude of the seismicity, and
M
min
=1
.
5
is the chosen cutoff magnitude. The estimated b-values in our dataset range from 0.6 to 1.2
along the central SAF.
Non-clustered seismicity rate and fraction of non-clustered seismicity:
Under the framework of Epidemic-type Aftershock-sequences Model (
51, 63, 64
), we as-
sume that the earthquake catalog is composed of two types of events in terms of temporal
clustering behavior: one is the clustered seismicity which is triggered by other events and their
occurrences follow the Omori law (
54, 64
). The other type is the non-clustered seismicity
which are independent from each other and follows the homogeneous Poisson process. The
former is equivalent to aftershocks, while the later represents the background seismicity. Our
analysis in the main text is based on the fraction of non-clustered seismicity in the earthquake
catalog. Based on the mild assumptions that general seismicity has Poissonian non-clustered
events superimposed with triggered aftershocks that follow the Omori law, it has been shown
that the non-clustered seismicity rate could be estimated with methods based on the distribution
of waiting times (interevent times),
t
, between consecutive earthquakes (
65
) proposed that the
probability density function of interevent times of general seismicity can be approximated as a
Gamma distribution:
P
(
t
)=
C
·
1
·
exp
i.e.,
t
(
)
where
C
is a normalization constant,
μ
is the non-clustered seismicity rate, and
is a model
parameter. The Gamma distribution,
(
)
, is assumed to be flexible enough to capture a range
of behavior that the interevent times of seismicity could exhibit. Molchan (
52
) and Hainzl et al.
(
41
), based on the theoretical analysis, further showed that the non-clustered seismicity rate,
μ
,
can be estimated from the first two sample moments of the interevent times:
μ
=
E
{
t
}
var
{
t
}
=
t
2
t
Here,
E
{
t
}
=
t
denotes the mean, and
var
{
t
}
=
2
t
is the variance of the interevent
times. Alternatively, we can have a similar form based on the normalized interevent times,
(where
=
t/
t
), relative to the mean value of a set of events (
41
):
p
(
)=
C
0
·
0
1
·
exp(
/
)
1
=
E
{
}
var
{
}
=
2
=
1
2
=
μ
where
is the total seismicity rate (
=1
/
t
). Therefore, the dimensionless factor,
1
/
encapsulates the fraction of non-clustered seismicity in the catalog, ranges from 0 to 1, gen-
erally. With
1
/
=1
, there are no triggered and temporally clustered earthquakes, and all
events are non-clustered events; While
1
/
=0
denotes all events are triggered and clustered
seismicity as a pure aftershock sequence that follows Omori law would exhibit. Intermediate
values in between represent the relative measure of non-clustered events activity in the catalog.
Note that the branching ratio (the fraction of triggered events in a catalog) is just the fraction of
non-clustered seismicity subtracted from one.
n
=1
1
2
[0
,
1]
Therefore, by calculating the mean of interevent times over the variance of interevent times,
we have an estimate of the non-clustered seismicity rate in the catalog. Further normalizing the
non-clustered seismicity rate with respect to the total seismicity rate gives us the fraction of non-
clustered events. Rather than relying on any particular triggering mechanism or declustering
algorithm, this technique is robust given the seismicity catalog is complete. It only depends on
a simple metric of interevent times from a set of events, and magnitude is not considered into
this formulation.
However, event subsets are not ideal under this framework if they deviate significantly from
the general mixture of independent non-clustered events (Poisson point process) and triggered
events (following the Omori law). For example, fitting a Gamma distribution to a perfect char-
acteristic repeating earthquake sequence with exactly the same interevent times would result
in a non-clustered rate of infinity, thus an infinite fraction of non-clustered events. Although
observed repeater sequences rarely recur with an exactly constant time interval, the high peri-
odicity of a repeater sequence can still result in a non-clustered event fraction larger than 1.0
based on this model, which is unphysical. A combination of many different repeating earth-
quake sequences (different interevent times) with many non-repeating events may be fine in
adopting this framework. Also, event subsets with only aftershocks are also shown to give an
overestimated non-clustered fraction (
32
).
Uncertainty of the estimates of the seismicity:
We use confidence intervals to present the uncertainty of estimates on seismicity properties.
The
95%
confidence interval of the non-clustered seismicity rate, the fraction of non-clustered
seismicity, and the b-value presented in Fig. S4, are based on the bootstrapping method. We
repeatedly and randomly resample the magnitude data for 10,000 times with replacement and
calculate the b-value each time. We also repeatedly resample the interevent-time data with
replacement and calculate the non-clustered seismicity rate and the fraction of non-clustered
seismicity. For these 10,000 sets of sampling distributions, we take the
2
.
5
th
and
97
.
5
th
per-
centiles to represent the lower bound and the upper bound of the
95%
confidence interval.
Supplementary Text
The non-clustered fractions in depth:
To better compare with fault slip models in two dimensions, we conduct the same analysis
in along-strike and depth cells (Fig. S5). We apply the calculation of non-clustered fractions to
cells on the fault plane with along-strike width of 2.5 km and a depth bin of 1.0 km. We omit
the cells with less than 20 events. The histogram counting the number of events in depth is also
binned into every 1 km.
The overall variation of the non-clustered fraction is consistent with the geodetic slip model
(
28
), where both ends are strongly coupled at nearly all depths, and the middle section is weakly
coupled. Limited by the seismicity distribution, we can only image the non-clustered fraction
at a narrow width at depth (mostly from 4 km to 12 km depth).
The consistency surprisingly exists in detail. The southern transition has a small patch with
a low non-clustered fraction at
195 km distance and 5–8 km depth (the cyan circle, A). This
patch corresponds to an asperity inferred from the geodetic model (
28
) 8 km below Parkfield. In
their model, a larger locked asperity has been inferred in the northern transition 40 km south of
San Juan Bautista. This asperity aligns with an area with low non-clustered fraction at 75–90 km
distance and 5–9 km depth (the purple circle, B). There is an asperity within the fast-creeping
segment at a distance
135 km and 7–10 km depth (red circle, C) caused by the 2012 M5.3
event at 8.6 km depth. This patch is consistent with the geodetic model as well (
27, 28
).
The depth-dependent non-clustered fractions are also consistent with the previously ob-
served microseismicity. The low non-clustered fractions within the north transition zone (the
purple circle) could be linked to the localized clusters (Fig. 2b and c in the main text) of in-
tense streaks of microearthquakes (
15, 28
). Yet, the paucity of events limits us from obtaining
estimates on the fault plane with even finer resolution and reasonable uncertainties.
The time-series of non-clustered fractions:
This section corresponds to Fig. S6 of this document. To ensure having enough events for
temporal analysis, we must sacrifice the spatial resolution with only four subsets of the entire
central SAF. The four subsets are: (1) the Loma Prieta segment; (2) the north transition zone;
(3) the fast-creeping segment; (4) the south transition zone near Parkfield. The blue curves in
Fig. S6 show the resulting time-series of the non-clustered fraction in four subsets, respectively.
The length of the moving time window is fixed at two years, with a step of two months.
As explained in the main text, there are issues when dividing the catalogs into fine time
windows due to a lack of events to fit a Gamma distribution of the interevent times properly.
Marsan et. al. (
32
) proposed a better time-series analysis, which we did not pursue in our study.
The basic idea is to adjust and extend the time windows adaptively to avoid only covering
aftershocks.
The effect of cutoff magnitude:
The estimated magnitude of completeness in our entire catalog is about M1.5 (Fig. S2–S3).
The interevent-time statistics below that cutoff magnitude is uncertain. Nevertheless, we can
explore the influences of choosing different cutoffs by incrementally increasing the cutoff mag-
nitude,
M
c
, and calculate the non-clustered fractions. The results are presented in Fig. S7
(
M
c
=2
.
0
) and Fig. S8 (
M
c
=2
.
5
). The general trend of non-clustered fractions along the
fault does not change. The fractions still correlate well with the surface creep rate regardless
of different cutoffs. But, by choosing a larger cutoff magnitude, the small aftershocks that are
excluded outnumber the mainshocks (non-clustered events). Using a larger cutoff magnitude
results in slightly higher values of non-clustered fractions, especially where small aftershocks
are populated. Nevertheless, the main point here is that the non-clustered fraction is indicative
to the coupling ratio. This conclusion is always obvious in the different magnitude ranges we
tested.
The unreleased moment within the north transition zone:
The transition zone north of the central creeping section may have the potential to produce
large earthquakes based on the intermediate level of coupling, which motivates an investigation
of the seismic moment release and accumulation. Assuming a conservative long-term plate
motion of 34 mm/yr (
29, 66
) and the averaged interseismic creep rate at 8 km depth of
9.8
mm/yr (
28
) along the transition zone, there is a slip deficit of 2.4 cm/yr. From 1984 to present
(37 years), this deficit corresponds to
3
.
0
10
19
N
·
m
of moment on a rectangular fault plane
which is 75 km long and 15 km wide, given a shear modulus of 30 GPa. Based on our dataset,
the released seismic moment by the seismicity within the transition zone adds up to only
3
.
3
10
17
N
·
m
, which is nearly two orders of magnitude smaller than the moment deficit. The
total unreleased moment for the north transition zone over 1984–2020 is equivalent to a
M
w
7.0
earthquake.
Fig. S1.
Hypocenters of the seismicity along the central SAF. The color is coded by event mag-
nitude. White stars mark the 1989 M6.9 Loma Prieta and the 2004 M6.0 Parkfield mainshock.
Red bars marked the 1989 Loma Prieta and the 2004 Parkfield rupture extent, respectively.
LP: Loma Prieta; SJB: San Juan Bautista; MR: Melendy Ranch; BW: Bitterwater; SC: Slack
Canyon; SAFOD: San Andreas Fault Observatory at Depth; Pk: Parkfield; Ch: Cholame.
Fig. S2.
Number of earthquakes each year from 1984 to 2021 (blue rectangles). The magnitude
of completeness estimated from the MAXC method are plotted with red dots. The values of
the magnitude of completeness range from 0.8–1.5. We thereby set a conservative minimum
magnitude
M
min
=1
.
5
for our entire analysis in this study.
Fig. S3.
Magnitude-frequency distribution of the earthquake catalog. Based on the assessment
on the magnitude of completeness, we specify a minimum magnitude,
M
min
=1
.
5
, which gives
us 24,082 events for our analysis over three decades from 1984 to 2021.
Non-clustered fraction
Non-clustered rate
Non-clustered fraction
Non-clustered rate
Fig. S4.
Non-clustered fraction, non-clustered rate, b-value, and the seismicity intensity along
the central SAF. Similar analysis to Fig. 2 of the main text. But the calculation here is only
applied to 14 separate along-strike windows, each 15-km wide, rather than a sliding window as
in Fig. 2. This is meant to show proper bootstrapping uncertainties with independent samples
between spatial windows. The red, blue, and orange shading shows the
95%
confidence interval
of the estimates. See Supplementary Text in the section ”Uncertainty of the estimates of the
seismicity”.
Fig. S5.
Non-clustered fraction estimated on a regular along-strike (2.5 km) and depth (1.0 km)
cells. The two white stars marked the hypocenters of the 1989 and the 2004 earthquakes. The
grey histogram on the right shows the depth distribution of events. Cells containing less than 20
events are colored as black to avoid biased interpretation. The colored circles indicate the fault
areas discussed in Supplementary Text in the section ”The non-clustered fractions in depth”.