## Abstract

Strain accumulated on the deep extension of some faults is episodically released during transient slow-slip events, which can subsequently load the shallow seismogenic region. At the San Andreas fault, the characteristics of slow-slip events are difficult to constrain geodetically due to their small deformation signal. Slow-slip events (SSEs) are often accompanied by coincident tremor bursts composed of many low-frequency earthquakes. Here, we probabilistically estimate the spatiotemporal clustering properties of low-frequency earthquakes detected along the central San Andreas fault. We find that tremor bursts follow a power-law spatial and temporal decay similar to earthquake aftershock sequences. The low-frequency earthquake clusters reveal that the underlying slow-slip events have two modes of rupture velocity. Compared to regular earthquakes, these slow-slip events have smaller stress drop and rupture velocity but follow similar magnitude-frequency, moment-area, and moment-duration scaling. Our results connect a broad spectrum of transient fault slip that spans several orders of magnitude in rupture velocity.

## INTRODUCTION

The establishment of continuous global positioning system (GPS) measurements led to the discovery of slow-slip events (SSEs) down dip of the seismogenically locked region of the Nankai (*1*) and Cascadia (*2*) subduction zones. Subsequently, tectonic tremors that correlate spatially and temporally with SSEs were found (*3*, *4*). These long-duration tremors have been inferred to be the superposition of many low-frequency earthquakes (LFEs) that represent asperities that are repeatedly driven to failure by surrounding aseismic slip (*5*). This interpretation that LFEs are markers of SSEs is supported by recent studies that managed to extract slow-slip deformation signal by using the timing of increased LFE rate as a guide to stack GPS time series (*6*, *7*).

Surface geodetic measurements can typically only detect SSEs above moment magnitude (*M*_{w}) ∼ 6 (*8*). Along the central San Andreas fault, tremors and LFEs have long been observed (*9*, *10*), but it was only recently that the deformation signature of *M*_{w} 4.9 SSEs were detected after stacking 20 such events using the timing of increased LFE rate (*7*). Since these SSEs in the deeper part of the fault might be episodically loading the shallow region that last ruptured in the 1857 magnitude 7.9 Fort Tejon earthquake (*11*), it is important that we have a more comprehensive catalog of these events. The bias toward only detecting the largest SSEs also limits our ability to robustly characterize their scaling properties. For instance, it is still debated whether the moment (*M*_{0}) of SSEs is proportional to their duration (*T*) (*12*) or follows the *T*^{3} scaling observed for regular earthquakes (*13*, *14*). A *M*_{0} ∝ *T* scaling relationship for SSEs could reflect fundamentally different underlying dynamics compared to regular earthquakes (*12*) or simply be the by-product of only cataloging the largest “bounded” events with aspect ratios >1 (*15*).

Previous studies have used tremors and LFEs to characterize SSEs invisible to geodetic measurements. Along the central San Andreas fault, Shelly (*10*) cataloged 88 LFE families (Fig. 1) that represent groups of LFEs with similar waveforms and, hence, similar source mechanisms and locations. Thomas *et al*. (*16*) clustered the LFEs of each individual family based on their recurrence intervals using empirically derived separation time scales, while Lengliné *et al.* (*17*) proposed a stochastic model that successfully reproduced the temporal clustering behavior of the LFEs. However, these studies did not account for interactions between LFE families and, hence, could not directly estimate the spatial properties of the underlying SSEs from the LFE clusters, even though the occurrence patterns of different families have been shown to correlate (Fig. 1C) (*18*, *19*). Clustering of tremors and LFEs that takes into account their spatial relationships has been attempted to identify secondary slip fronts at subduction zones but relied on rather ad hoc definitions of what constitutes a cluster (*20*, *21*).

Here, we estimate the spatiotemporal clustering properties of LFEs detected along the central San Andreas fault (*10*) probabilistically with minimal model assumptions and a priori parameterization. We then use the extracted LFE clusters to estimate the area, average slip, moment, duration, stress drop, and rupture velocity of the underlying SSEs and explore their scaling properties.

## RESULTS

### Clustering properties of LFEs

We analyze the catalog of 88 LFE families that span ∼150 km along the central San Andreas fault (Fig. 1) (*10*), which includes more than 1 million events from 2001 to 2016. We limit our analysis to between 2006 and 2016 to minimize the impact of the 2004 *M*_{w} 6 Parkfield earthquake. This leaves us with a catalog of ∼750,000 LFEs. The most populous family has ∼22,000 events, while the least populous family has ∼2500 events over this time period. We model the LFE rate at time *t* as*D* is the number of LFE families, μ* ^{x}* is the uniform background rate of family

*x*,

*K*encodes the excitation strength, i.e., the number of events in family

^{xy}*x*“excited” by an event in family

*y*on average,

*g*is the normalized time-dependent excitation kernel, and

*j*from family

*y*. Therefore, we are only assuming that (i) the LFE clustering behavior is linear, i.e., the contribution of different LFEs can be added up; (ii) a mean-field response to the occurrence of an LFE, i.e., two LFEs from the same family are modeled similarly; and (iii) the excitation kernel

*g*(

*t*) does not vary between different families, although we obtain similar results when relaxing this assumption (see Materials and Methods). We do not make an a priori assumption about the shape of

*g*(

*t*) but instead discretize it as piece-wise constant. We then adopt an expectation-maximization approach (

*22*,

*23*) to estimate the parameters μ

*,*

^{x}*K*, and

^{xy}*g*(

*t*) (see Materials and Methods). We verified the algorithm’s ability to correctly estimate the model parameters using a synthetic catalog (see Materials and Methods).

*g*(*t*) characterizes the temporal clustering property of the LFEs. We find that *g* follows a power-law decay with time up to 10 days, with a plateau between 0.02 and 0.2 days (Fig. 2A). We fit the decay in the range of [2 × 10^{−4} to 2 × 10^{−2}] days and between 0.2 and 10 days separately and obtain *g* ∝ *t*^{−1.8} for both time ranges. The shape of *g*(*t*) and the power-law exponent that we obtain are similar to what Lengliné *et al.* (*17*) obtained from applying a comparable stochastic model to one LFE family at a time, as well as the stacked interevent time density of the different LFE families (fig. S1). Therefore, the power-law decay of *g*, which is similar to the temporal evolution of earthquake aftershock sequences often referred to as the Omori-Utsu law (*24*), is a robust feature of the LFE catalog. Lengliné *et al.* (*17*) concluded that this feature is unlikely to have arose from direct triggering of LFEs by the stress change due to a preceding LFE since the number of excited LFEs is not correlated with the amplitude of “mainshock” LFEs in their analysis. Therefore, the power-law decay of *g* likely reflects the triggering of LFEs by changes in stress or loading rate due to underlying SSE. The plateau of *g* between 0.02 and 0.2 days (30 min and 5 hours) is rather perplexing but indicates that there are two clustering time scales. When we allow *g*(*t*) to vary between LFE families (see Materials and Methods and fig. S2A), we find that *g*(*t*) has such a plateau only for LFE families that have short-duration bursts occurring within long-duration bursts, i.e., a trimodal interevent time distribution (fig. S3) (*16*), similar to secondary slip fronts observed within SSEs at the Cascadia and Nankai subduction zones (*21*, *25*).

*K ^{xy}* characterizes the excitation strength between LFE families and thus the spatial clustering property of the LFEs. We find that a family that was previously suggested to be isolated (Fig. 1) (

*26*) has the second smallest interfamily

*K*value, lending confidence that

*K*indeed captures the interaction between LFE families. On average,

*K*decays with interfamily distance (Fig. 2B), even though interfamily distance is not present in our model (Eq. 1). The decay of

*K*with distance is consistent with previous observations that the occurrence pattern of nearby LFE families is correlated (

*18*,

*19*). For along-strike distance, the decay of

*K*is well approximated by a power law up to 16 km, above which the value of

*K*saturates possibly due to the values being too small to be resolvable (Fig. 2B). We fit the power-law decay between 1 and 16 km and obtain

*K*∝

*d*

^{−2.8}. For along-dip distance, we obtain

*K*∝

*d*

^{−2.5}. The fast decay of

*K*with distance is consistent with Trugman

*et al*. (

*19*) obtaining groups that typically span <15 km when grouping LFE families based on the similarity of their long–time scale occurrence pattern. The power-law exponent of ∼3 is similar to the expected earthquake aftershock density decay with distance from the mainshock if aftershocks are triggered by static stress change. However, a model where LFEs are triggered by static stress change due to SSEs is inconsistent with interfamily excitation in the along-strike direction being 10 times stronger compared to the along-dip direction (Fig. 2B). Therefore, the decay of

*K*with distance probably reflects the likelihood of the underlying SSEs to grow spatially, with there being a greater tendency to propagate along strike than along dip. This favors the model where LFEs are small asperities driven to failure by surrounding aseismic slip. Lastly, we find that the decay of

*K*is similar for both southeast (the excited family is located southeast of the exciting family) and northwest excitation (Fig. 2B). This suggests that the underlying SSEs are equally likely to propagate in either direction, even though previous studies have suggested that earthquakes along the Parkfield segment of the San Andreas fault preferentially rupture to the southeast [e.g., (

*27*)].

### Estimating properties of underlying SSEs

Now that we have estimated μ, *K*, and *g*, which govern the LFE rate at any given time, for each LFE *i*, we can calculate the probability that it is a background event and the probability that it was excited by each preceding LFE *j* (see Materials and Methods). Using the probabilities associated with each LFE, we perform stochastic clustering (*28*) of the events to isolate individual LFE bursts. Each cluster includes one background event and the events it directly and indirectly “excites.” We interpret this background event to mark the initiation of an SSE, while the excited events reflect the subsequent evolution of the SSE. Stochastic clustering is not unique, and each catalog that it produces represents a sampling of the underlying structure. Here, we present a catalog from one iteration, which includes 16,327 LFE clusters that involved ≥2 LFE families, but have verified that the presented statistics remain stable between different iterations as expected. For each of these clusters, we calculate their along-strike extent (*L*) and depth extent (*W*). We then infer the underlying SSE to have a rupture area *A* = *LW* (fig. S4). The time difference between the first and last LFEs in a cluster is taken as the SSE duration (*T*). We then estimate the average rupture velocity *5*) and our observed clustering properties of LFEs appear to support this model, our estimating these SSE properties using LFE clusters only requires that they are spatially and temporally coincident (*7*) and is valid regardless of the physical mechanism behind their coincidence.

Thomas *et al.* (*16*) have shown that LFEs on the San Andreas fault can be used to meter slip. The long-term slip rate of the 150-km-long creeping segment of the San Andreas fault is approximately 34 mm/year [e.g., (*29*)]. For each LFE family *x*, we estimate the slip per LFE (*d _{x}*) by multiplying the long-term slip rate by the catalog duration (10 years) before dividing by the total number of LFEs in that family. For each LFE cluster, we can then estimate the total slip of each LFE family

*x*by multiplying

*d*by the number of LFEs of that family that is part of the cluster. We assume the total slip of each LFE family represents a point sample of the slip distribution over the rupture area of the underlying SSE, i.e., there is slip throughout the area delineated by each LFE cluster, but we only have an estimate of the slip amount where the LFE families are located. For each LFE cluster that involved ≥2 LFE families, we then take the mean slip of all the activated LFE families as the average slip

_{x}*M*

_{o}estimates by a constant factor and would not affect its inferred scaling relationship with other SSE properties. We calculate

*M*

_{w}of 5.1. This is just slightly larger than the geodetically determined average

*M*

_{w}of 4.9 for the SSEs associated with the 20 largest bursts of 10 LFE families on the northwest Parkfield segment of the San Andreas fault (

*7*), suggesting that our derived SSE moments are reasonable.

## SSE scaling properties

Earthquakes follow a power-law magnitude-frequency distribution given as log_{10}(*N*) = *a* − *bM*, where *N* is the number of earthquakes of magnitude ≥ *M*, and *a* and *b* are constants (*30*). We find that the SSEs also approximately follow a power-law magnitude-frequency distribution (Fig. 3A). We obtain a *b* value of 1.61 ± 0.04 (*31*, *32*) when taking a magnitude of completeness *M*_{c} = 3.9. Nonvolcanic tremors that are approximately spatially coincident with our LFE families in the southeast have been suggested to also follow a power-law magnitude-frequency distribution with a *b* value of 2.50 (*33*). The same study found that regular earthquakes on the creeping section of the San Andreas fault, which are on the shallower portion of the fault above our LFE families in the northwest, have 0.8 < *b* < 2.0. SSEs at the Cascadia subduction zone have also been suggested to follow a power-law magnitude-frequency distribution with *b* values between 0.8 and 1.0 (*14*, *34*). The differences in *b* values might reflect differences in stress condition (*35*, *36*) and fault roughness (*37*).

We find that the SSE moment and rupture area approximately follow the *M*_{o} ∝ *A*^{1.5} scaling of regular earthquakes (Fig. 3B). Since we estimated *M*_{o} from *A* and *M*_{o} < 10^{13.5} (Fig. 3B), probably due to the poorer rupture area estimates for these smaller events since their length and width are comparable to the LFE location uncertainties of 1 to 2 km (*10*). For bins of *M*_{o} > 10^{13.5} N · m, we calculate the median rupture area. Fitting these median values, we obtain *M*_{o} ∝ *A*^{1.5} (Fig. 3B). This is consistent with our SSEs being mostly “unbounded” events (*15*) where both the length and width are still growing as the rupture propagates (fig. S5B), although these events show a greater tendency to rupture in the along-strike direction (Fig. 2B) with the median SSE aspect ratio *W* eventually saturates while *L* and *W* of ∼13 km (fig. S5B) could reflect the maximum channel width of SSEs at the central San Andreas fault, as determined by rheological change with depth, or simply a coincidence of asperities’ distribution that resulted in the LFEs being located in a narrow 13-km band (Fig. 1B) even though SSEs can rupture beyond this zone.

The SSE stress drops are mostly within the range of 1 to 10 kPa (Fig. 3B), with a median value of 6 kPa based on a circular crack model (*38*). This is a few orders of magnitude smaller than the stress drop of regular earthquakes on the San Andreas fault, which are estimated to be on the order of 0.1 to 100 MPa (*39*). The SSEs have a median stress drop of 2 kPa when assuming strike-slip faulting on a rectangular fault instead (*40*). Our results are comparable to the estimated stress drop of ∼10 kPa from spectral analysis of two LFE families on the San Andreas fault (*41*) and is consistent with the low stress drops (0.1 to 10 kPa) typically obtained for slow earthquake phenomena [e.g., (*12*, *14*, *21*, *34*, *42*)].

The logarithm of the SSE duration has a bimodal distribution (Fig. 3C) with a clear separation at approximately 10^{3.5} s (84 min), which is in the time range when *g* has a plateau (Fig. 2A). This translates to a bimodal distribution in the rupture velocity, with two modes at 6 and 708 km/day (fig. S6). These velocities are within the range of tremor migration velocities previously observed at the San Andreas fault (*18*), as well as at the Nankai (*5*) and Cascadia subduction zones (*21*). We use *T* = 10^{3.5} s as a boundary to separate the SSEs into two populations. For the shorter-duration population, we calculate the median *T* for *M*_{o} bins in the range of [10^{11} to 10^{15}] N · m. For the longer-duration population, we calculate the median *T* for *M*_{o} bins in the range of [10^{12.5} to 10^{16.5}] N · m. Fitting these median values, we obtain *M*_{o} ∝ *T*^{3.1} and *M*_{o} ∝ *T*^{2.8} for the shorter- and longer-duration SSE populations, respectively (Fig. 3D). Recently, Michel *et al*. (*14*) suggested that 40 *M*_{w} 5.3 to 6.8 SSEs cataloged by inverting GPS measurements at the Cascadia subduction zone likely follow a *M*_{o} ∝ *T*^{3} relationship. However, their best-fit value is *M*_{o} ∝ *T*^{5}, and for these large SSEs that rupture the entire down-dip width of the slip zone, a linear scaling is expected instead (*15*). Therefore, questions remain regarding how well constrained the scaling relationship is given the small sample size, as well as whether the observed scaling relationship is unique to the Cascadia subduction zone. At the Mexican subduction zone, Frank and Brodsky (*13*) calibrated daily median *S* wave amplitudes of LFEs from 58 families to GPS-measured moment rate of 12 SSEs. They then use the calibrated seismic to geodetic moment rates scaling to estimate daily moment release from the LFE catalog and find that the underlying SSEs follow a *M*_{o} ∝ *T*^{3} relationship. Since we obtain the same scaling relationship for SSEs at a transform fault using a fundamentally different method, it is likely a universal property that unbounded SSEs follow a *M*_{o} ∝ *T*^{3} scaling (*15*) similar to regular earthquakes (*38*). We also show that there are two populations of SSEs at the San Andreas fault with vastly different rupture velocities, with each population approximately following a cubic moment-duration scaling (Figs. 3D and 4). Therefore, the previously apparent *M*_{o} ∝ *T* scaling suggested to be a unique characteristic of slow earthquakes (*12*) is possibly an artifact of fitting the relationship across different populations of slow-slip phenomena spanning a wide range of rupture velocities.

## CONCLUSIONS

We have shown that LFEs can be used to catalog the large number of SSEs that are episodically loading the shallow seismogenic zone. We find that the LFEs’ clustering properties are similar to earthquake aftershock sequences but with interfamily excitation in the along-strike direction being 10 times stronger than in the along-dip direction. The underlying SSEs have two modes of rupture velocity with stress drop and rupture velocity a few orders of magnitude smaller than regular earthquakes. However, the SSEs follow similar magnitude-frequency, moment-area, and moment-duration scaling as regular earthquakes, suggesting that transient fault slip with velocities spanning many orders of magnitude may be governed by the same dynamics. Our observations provide important constraints on the relationships between LFEs, SSEs, and regular earthquakes.

## MATERIALS AND METHODS

### Modeling LFE rate

For the most generalized model, the LFE rate at time *t* can be modeled as*D* is the number of LFE families, μ* ^{x}* is the uniform background rate of LFE family

*x*,

*g*is the time-dependent excitation kernel that encodes the influence of an LFE from family

^{xy}*y*on the future rate of family

*x*, and

*j*from family

*y*.

*g*(

*t*) can be discretized as piece-wise constant

*T*are the discretization times and

_{m}*m*∈ [1 :

*M*] with

*T*

_{1}= 0 and

*T*= 10 days. For

_{M}*M*= 20 and

*D*= 88, we would have to solve for ∼155,000 parameters. Therefore, we simplify the model by assuming that

*g*(

*t*) does not vary between different families such that

*K*encodes the excitation strength, i.e., the number of events in family

^{xy}*x*excited by an event in family

*y*on average, and

*g*(

*t*) is normalized such that it represents a probability density function

*t*is the discretization time step of

_{m}*g*, i.e., δ

_{m}*t*=

_{m}*T*

_{m + 1}−

*T*. With this model, we only have to solve for ∼8000 parameters. ∑

_{m}_{x, y}

*K*has to be <

^{xy}*D*for the model to be stable; ∑

_{x, y}

*K*>

^{xy}*D*would result in an infinite number of events within a finite time period, i.e., λ growing to infinity. We adopt an expectation-maximization approach to estimate μ

*,*

^{x}*K*, and

^{xy}*g*(

_{m}*22*,

*23*). The expectation step involves computing the probability that event

*i*from family

*x*was excited by event

*j*from family

*y*

*t*, and the probability that event

_{i}*i*from family

*x*was a background event (not excited by any previous events)

The log-likelihood function associated with Eq. 4 is thus*T* is the duration of the time series, *n _{y}* is the number of events from family

*y*.

The maximization step then involves updating the background rates as*g* following Eq. 5. We start with initial guesses of μ* ^{x}* (uniform values of 1 for all families),

*g*(random values between 0 and 1), and

_{m}*K*(random values between 0 and 1). We then iterate through the expectation-maximization steps until the difference of the estimated parameters between two successive iterations is smaller than a certain threshold. We obtain ∑

^{xy}_{x, y}

*K*= 81.6. We generate a 10-year synthetic catalog with the

^{xy}*g*, μ, and

*K*that we obtained for the San Andreas LFE catalog, with

^{xy}*K*first multiplied by 0.8 to test whether we can resolve parameters from a catalog with weaker interevent excitation. We find that we can correctly estimate the model parameters from the synthetic catalog using the proposed expectation-maximization algorithm (fig. S7).

^{xy}Since it has been observed that different LFE families have different interevent time distribution (*16*), we explored having one *g*(*t*) per family such that*g ^{x}* is the time-dependent excitation kernel that encodes the influence of an LFE on the future rate of family

*x*. With this model, we have to solve for ∼9500 parameters and obtain ∑

_{x, y}

*K*= 83.3. While

^{xy}*g*(

^{x}*t*)’s shape varies depending on the LFE family’s interevent time distribution (fig. S3), we obtain similar scaling properties (figs. S2B and S8) as when using Eq. 4, except that we now obtain a

*b*value of 1.35 (fig. S8A) and

*M*

_{o}∝

*T*

^{4.2}for the shorter-duration SSE population (fig. S8D).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/33/eabb2489/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank B. Yuan, C. Cattania, D. M. Saffer, G. C. Beroza, P. Segall, W. L. Ellsworth, and W. B. Frank for fruitful discussions. We also thank B. Rousset, C. H. Scholz, and H. Houston for providing helpful feedback on the manuscript.

**Funding:**Y.J.T. acknowledges support by the Chateaubriand fellowship and the “Make Our Planet Great Again” initiative.

**Author contributions:**Y.J.T. and D.M. designed the study. Y.J.T. performed the data analysis, made the figures, and wrote the initial draft of the manuscript. All authors contributed to the interpretation of the results and writing of the manuscript.

**Competing interest:**The authors declare that they have no competing interests.

**Data and materials availability:**The LFE catalog is available from the supplementary materials of Shelly (

*10*). The statistical modeling was performed using a modified version of an open-source software (

*43*) that can be downloaded from https://x-datainitiative.github.io/tick/. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper can be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).