Dans les deux études suivantes, nous avons chercher à observer les détails de la dynamique d'un séisme lent dans la région des Cascades pour répondre à ces questions.

## 2 Les fronts secondaires de glissement associés aux séismes lents

Le glissement lent sur la faille est généralement imagé par inversion de données GNSS. La résolution de ces données et les limites des techniques d'inversion conduisent à imager des distributions de glissement très lisses spatialement et temporellement. La micro-sismicité qui leur est associée présente, au contraire, une très bonne résolution temporelle et une assez bonne résolution spatiale. Si l'ont fait l'hypothèse que tremors et LFEs sont les manifestations passives de la rupture de micro-aspérités déstabilisées par le glissement lent sous-jacent [Shelly et al., 2007], on peut considérer que là où l'on observe tremors et LFEs, il y a glissement lent, et ainsi avoir indirectement accès à une vision du glissement beaucoup mieux résolue que par inversion de données GNSS. Pour analyser la dynamique des séismes lents, nous avons ainsi, dans l'étude qui suit, identifié, de manière automatique, un grand nombre de migrations de tremors/LFEs le long de la faille de subduction des Cascades. Nous avons interprété ces migrations d'évènements micro-sismiques comme les témoins de fronts de glissement lent secondaires – que nous avons nommé *Secondary Slip Fronts* (SSFs) – au sein des séismes lents. Nous avons ainsi mis en évidence que les séismes lents ont une dynamique spatio-temporelle très complexe, composée de nombreux sous-évènements. Ces sous-évènements, les SSFs, semblent par ailleurs – même si les approximations faites pour estimer leurs magnitudes sont susceptibles de biaiser les estimations vers ce résultat – obéir à la loi d'échelle proposée par *Ide et al.* [2007] (magnitude proportionnelle à la durée) pour les SSEs, les LFEs et les VLFEs (Figure I.6) dans une gamme de magnitude et de durée inaccessible à la géodésie et à la sismologie conventionnelles.

Cette étude a été réalisée à l'Université d'Oregon, en collaboration avec Amanda M. Thomas, Jessica C. Hawthorne, Robert M. Skarbek, Alan W. Rempel and Randy D. Krogstad. Elle a été publiée en 2017 dans le journal *Earth and Planetary Science Letters* [Bletery et al., 2017b].
<figure><img src="image_35.png" /></figure>

ELSEVIER

Contents lists available at ScienceDirect

Earth and Planetary Science Letters

www.elsevier.com/locate/eps1

<figure><img src="image_36.png" /></figure>

# Characteristics of secondary slip fronts associated with slow earthquakes in Cascadia

<figure><img src="image_37.png" /></figure>

CrossMark

Quentin Bletery<sup>a,*</sup>, Amanda M. Thomas<sup>a</sup>, Jessica C. Hawthorne<sup>b</sup>, Robert M. Skarbek<sup>c</sup>,
Alan W. Rempel<sup>a</sup>, Randy D. Krogstad<sup>a</sup>

## ARTICLE INFO

**Article history:**

Received 10 November 2016

Received in revised form 30 January 2017

Accepted 31 January 2017

Editor: P. Shearer

**Keywords:**

slow slip

low-frequency earthquakes

nonvolcanic tremor

secondary slip fronts

Cascadia

## ABSTRACT

We implement an algorithm to automatically detect migrations of low frequency earthquakes at time scales between 30 min and 32 h during the 2003, 2004 and 2005 slow slip events in Cascadia. We interpret these migrations of seismicity as a passive manifestation of secondary slip fronts (SSFs) that propagate faster than the main front. We identify the dominant features of 383 SSFs, including time, location, duration, area, propagation velocity and estimate: their moment, stress drop, slip, and slip rate. We apply the same algorithm to continuous tremor detection in Cascadia between 2009 and 2015 and characterize 693 SSFs at time scales between 4 h and 32 h. We identify — to our knowledge for the first time — numerous 11–22.5 h long SSFs that propagate at velocities intermediate between slow slip events and previously reported SSFs. The systematic detection of SSFs fills a gap between seismically and geodetically detectable slow earthquake processes. Analysis of SSF basic features indicates a wide range of stress drops and slip rates (with medians of 5.8 kPa and 1.1 mm/h) as well as an intriguing relationship between SSF direction and duration that was observed in other contexts and could potentially help discriminate between the different physical models proposed to explain slow slip phenomena.

© 2017 Elsevier B.V. All rights reserved.

## 1. Introduction

Slow slip events (SSEs) in subduction zones propagate along the plate interface at velocities of 5–10 km/day and are largely confined to the region known as the transition zone, located down-dip of the seismogenically locked zone (Dragert et al., 2001; Bartlow et al., 2011). As SSEs propagate, small on-fault asperities capable of generating seismic radiation fail in earthquake-like events known as low-frequency earthquakes (LFEs) and nonvolcanic tremor (Obara, 2002; Kao et al., 2006; Ide, 2010; Ueno et al., 2010; Houston et al., 2011). Multiple lines of evidence indicate that LFEs and tremor — which has been interpreted as a superposition of LFEs (Shelly et al., 2007b) — are a passive manifestation of surrounding aseismic slow slip (Shelly et al., 2007b). The occurrence of LFEs and tremor during SSEs are not random. Instead, spatially and temporally coherent migrations of LFEs and tremor can often be identified within the actively slipping portion of the fault after the main slip front associated with the SSE has passed. These

migration patterns are thought to represent secondary slip fronts (SSFs) occurring over several different lengths and timescales. For example, Houston et al. (2011) and Royer et al. (2015) used tremor and LFEs to image approximately 40 km long slip fronts — known as rapid tremor reversals (RTRs) — that propagate back in the opposite direction but 10 to 40 times faster than the main slip front. Ghosh et al. (2010) used tremor in Cascadia to image seismic migrations — known as streaks — that propagate back and forth along the main slip front (i.e., perpendicular to its propagation direction) but approximately 100 times faster. Rubin and Armbruster (2013) reported tremor migrations at the kilometer scale that propagate along strike either towards or away from the main front, but approximately ten times faster. Finally, Peng et al. (2015a) and Peng and Rubin (2016) found evidence for complex tremor migration patterns that evolve as the main front propagates. While all of the aforementioned studies were conducted in Cascadia, similar observations have been reported in Japan (Shelly et al., 2007a; Ueno et al., 2010; Obara et al., 2012), Mexico (Radiguet et al., 2012; Frank et al., 2014; Peng and Rubin, 2017) and along the San Andreas fault (Ryberg et al., 2010; Shelly, 2009, 2010; Peng et al., 2015b), suggesting that SSFs may be a ubiquitous feature of slow slip.

<footnote><sup>a</sup> Department of Earth Sciences, University of Oregon, 1272 University of Oregon, Eugene, OR 97403, USA</footnote>

<footnote><sup>b</sup> University of Leeds, School of Earth and Environment, Leeds, United Kingdom</footnote>

<footnote><sup>c</sup> Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY, USA</footnote>

<footnote>* Corresponding author.
E-mail address: qbletery@uoregon.edu (Q. Bletery).</footnote>
<figure><img src="image_38.png" /><figcaption>**Fig. 1.** Algorithm scheme. See main text for description.</figcaption></figure>

Slow-slip phenomena require the fault slip rate to increase to observed speeds, typically 10 to 100 times the plate rate, but abstain from accelerating to typical earthquake slip rates. The specific physical mechanisms responsible for imposing this speed limit are still debated. SSFs generally propagate at least an order of magnitude faster than the main front associated with the SSE, and estimates of the ratio of the SSF stress drop to that of the main front range between one tenth and one (Rubin and Armbruster, 2013; Royer et al., 2015; Hawthorne et al., 2016). These stress drop estimates combined with the propagation speeds mentioned above imply that SSFs must be slipping either as fast as or faster than the main front. Hence, reproducing stress drops, slip and propagation speeds of both the main and secondary fronts requires a rheology that limits the speed of the main front while allowing the SSFs to slip one to two orders of magnitude faster. Such behavior is difficult to reproduce with a single state variable friction law (Rubin, 2011). Additionally, SSFs appear to be sensitive to very small (≈ 1 kPa) tidal stress perturbations (Thomas et al., 2013; Royer et al., 2015; Houston, 2015) allowing them to be used to probe in situ conditions in deep faults. Observations of SSFs, in addition to SSEs, may provide useful constraints on physical models of slow slip processes (Rubin, 2011; Hawthorne and Rubin, 2013).

To date, most studies that have documented SSFs have relied on subjective identification methods such as visual identification, often using only one dimension of space and time (Shelly et al., 2007a; Ueno et al., 2010; Ghosh et al., 2010; Houston et al., 2011; Rubin and Armbruster, 2013; Royer et al., 2015; Peng et al., 2015a; Savard and Bostock, 2015). While these approaches successfully find some SSFs, such subjective identification makes it possible and even likely that many SSFs go undetected. Here we develop an algorithm that systematically identifies SSFs based on spatial and temporal migration of tremor/LFEs and discuss its advantages compared to the similar method developed by Obara et al. (2012). We apply our algorithm to both tremor and LFE datasets resulting from SSEs that occurred between 2003 and 2015 in Cascadia. We detect over 1000 SSFs, determine their locations, duration, directions, and estimate some of their physical properties, such as moment, average slip, stress drop and slip rate.

## 2. Data and algorithm

We use two types of data: a catalog of LFEs located near Vancouver Island during the 2003, 2004 and 2005 SSEs (Savard and Bostock, 2015) and a catalog of continuous tremor detections from August 2009 to December 2015 that occur throughout the Cascadia subduction zone (Wech and Creager, 2008; Wech, 2010). The LFE catalog was obtained using cross-station waveform correlation; thanks to a dense campaign seismic network, it has very good spatial (location errors of < 1 km) and temporal resolution (Savard and Bostock, 2015). The tremor catalog contains locations and detection times of tremor events determined using waveform envelope correlation and clustering (Wech and Creager, 2008; Wech, 2010). Detections correspond to 5-minute seismogram windows that likely represent a superposition of LFEs (Shelly et al., 2007b). Therefore, resolution of both time (each detection is separated by at least 2.5 min) and location (5–10 km) is worse than for the LFE catalog. Because the tremor detection algorithm is independently implemented on several adjacent seismic networks, some tremor detections appear twice in the final catalog (Wech, 2010). To eliminate duplicate detections, when two tremors are detected at the same time within 25 km, we exclude one (Wech, 2010).

These two catalogs provide a large number of identified LFE and tremor events. To identify possible SSFs, we analyze LFE and tremor migration patterns in 3 dimensions: distance along strike, distance along dip, and time. The algorithm is illustrated in Fig. 1. As in the approach of Obara et al. (2012), we look for SSFs in time windows of various durations. The different time windows considered in our analysis are 30 min, 1 h, 2 h, 4 h, 8 h, 16 h, 32 h for the LFE catalog and 4 h, 8 h, 16 h, 32 h for the tremor catalog. As individual detections in the tremor catalog are all separated by at least 2.5 min (Wech and Creager, 2008), the number of potential detections in a $\le$ 2 h time window is too small to analyze their potential migration. The analysis is performed independently on the different time windows. Therefore, tremor/LFE detections may theoretically be part of more than one SSF in different time windows; thus a short duration SSF may be a part of a longer duration SSF.

We first identify clusters of LFE and tremor detections in the individual catalogs. We use the subclust Matlab routine to identify
LFE/tremor clusters of different durations. The subclust routine estimates the centers of subclusters of a given duration in a data set using subtractive clustering (Chiu, 1994; Yager and Filev, 1994). We set the duration of the subclusters to correspond to the different time windows. The subclust routine then returns the tremor/LFE clusters of 30 min to 32 h duration with the highest concentration of events. For each identified cluster, we then calculate the linear regression of the distance along strike as a function of time and exclude points farther than two standard deviations ($2\sigma$) from the regression line. After removing those points, we do the same linear regression but for the distance along dip as a function of time and remove events farther than $2\sigma$ from the linear regression. Finally, we project the remaining points along the best-fitting migration axis and remove again events that lie outside $2\sigma$ from the linear regression of the distance along the migration axis as a function of time. The remaining points are then tested: if the cluster still contains more than 20 LFE or tremor detections and if the root mean square (RMS) of their distance to the last regression is smaller than 15% of the migration length, then the cluster is identified as an SSF. For a given time window, if the SSF length is small, then the front is not moving quickly (the signal is low). Thus we require the deviation from the propagation axis (the noise) to be very low (to ensure a decent signal to noise ratio). If the SSF length is larger, then the front is clearly moving (the signal is larger) and we can be more tolerant with the deviation from the propagation axis (as we have a clear signal). We repeat the process for all the clusters. We then remove the LFE or tremor detections identified as part of an SSF from the catalog and rerun the clustering routine on the new catalog with identified SSFs removed. We iterate until the catalog no longer contains any LFE or tremor. If no SSF is identified on one iteration, we remove all the tested clusters from the catalog and iterate once more in the truncated catalog. All the detected SSFs are then tested to ensure that each SSF is continuous rather than composed of two distinct SSFs. We subdivide the successfully identified SSFs into 4 equal portions along their best-fitting migration axis and 3 equal periods and if one of the parts or periods contains less than 1/20th of the number of LFE/tremor detections that compose the SSF, the SSF is rejected.

Our approach is inspired by Obara et al. (2012) (who implemented an automatic SSF detection based on principal component analysis of tremor in Japan), but has some significant differences. Obara et al. (2012) identified SSFs using moving time windows shifted by discrete increments. As it is unlikely that SSFs all start and end at the specific times of arbitrary time windows, this may result in detecting SSFs that either are just a part of a larger undetected SSFs or include LFE/tremor that should not be associated with the SSF, potentially affecting the apparent migration velocities. The use of a clustering subroutine allows for flexibility in the start and end time of the explored SSFs — as tested time windows are not shifted by discrete increments but determined in respect of high concentrations of tremor/LFEs — substantially increasing the number of SSF detections. Another significant difference is the detection thresholds. In the work of Obara et al. (2012), SSF acceptance is based on linearity, angular difference and in the number of tremor detections, all of which vary with the given time window. We apply only two thresholds, which are the same for all time windows: the number of LFE/tremor detections (20 for all time windows) and a deviation criteria based on the ratio of the misfit to the SSF length (required to be ≤ 15%), which can be interpreted as a signal to noise ratio. Our choices of 20 LFE/tremor detections and 15% are arbitrary — as are the thresholds used by Obara et al. (2012) — but these values do provide visually coherent patterns in all the detected SSFs (Figs. 2, 3).

To ensure that our algorithm does not interpret noise as SSF, we test it on synthetic catalogs. To produce the synthetic catalogs,

we use the actual catalog locations and randomize the time of occurrence. We run the algorithm for all the time windows of our analysis and detect no SSF in the synthetic catalogs.

### 3. Spatial distribution of secondary slip fronts

We identified 383 SSFs in the LFE catalog of the 2003, 2004, and 2005 SSEs beneath Vancouver Island and 693 SSFs in the tremor catalog between August 2009 and December 2015 in Cascadia. While SSFs lasting from 10 min to 11 h have been documented in multiple subduction zones (Ghosh et al., 2010; Houston et al., 2011; Obara et al., 2012; Rubin and Armbruster, 2013; Royer et al., 2015; Peng et al., 2015a; Savard and Bostock, 2015), we identified — to our knowledge for the first time — numerous SSFs lasting between 11 and 22.5 h in both catalogs (243 in the tremor catalog, 57 in LFEs catalog). They propagate at intermediate speeds (≈ 3 km/h) between those of the main front (≈ 0.2–0.5 km/h) and previously documented SSFs (≈ 40 km/h) (Ghosh et al., 2010; Houston et al., 2011; Obara et al., 2012; Rubin and Armbruster, 2013; Royer et al., 2015; Peng et al., 2015a; Savard and Bostock, 2015).

The spatial distribution of SSFs is of particular interest because it can reveal atypical behaviors of some fault areas and may highlight spatial variations of fault properties. Fig. 4 shows the distribution of tremor derived SSFs during the 2009–2015 period along the Cascadia megathrust by time window. Note that some long duration tremor derived SSFs may correspond to the main front of short duration SSEs. In the inset diagrams we plot the number of SSFs as a function of their propagation directions. Some areas show a clear tendency to host numerous SSFs (for instance at latitude 40°N–41.5°N, 42°N–42.5°N, or 44.5°N–45.5°N) whereas other regions had no detectable SSFs. Additionally, there is a clear preference for along-strike propagation for the majority of the tremor derived SSFs. Both the spatial distribution and propagation directions are consistent from SSE to SSE, suggesting that the physical mechanism responsible for generating SSFs does not change over a timescale of several years.

Because both the temporal and spatial resolution are much better in the LFE catalog, we are able to detect short timescale SSFs beneath Vancouver island (Fig. 5). Since the 3 SSEs that drove LFEs in this area all roughly propagated in the same direction (blue arrow in Fig. 5), this catalog is also more convenient to analyze the direction of SSFs with respect to the direction of SSEs. SSFs propagate in all directions and encompass phenomena previously described as RTRs or streaks. They also encompass slip fronts that propagate in the direction of the SSE (Fig. 5). Such forward-propagating SSFs are hard to distinguish by eye from LFE or tremor migration associated with the main front. They typically initiate after the main front has passed and migrate towards the tip of the main front. Interestingly, we observe a clear relationship between the direction and the duration of SSFs. The vast majority of short duration (i.e. 30 min and 1 h) SSFs propagate along dip. As the SSF duration increases, the proportion of SSFs propagating along strike gradually increases (see Rose diagrams in Fig. 5) resulting in the longest duration LFE derived SSFs propagating almost solely along strike. Similar relationships were observed by Ghosh et al. (2010) and Houston et al. (2011) in Cascadia and by Obara et al. (2012) beneath the Kii peninsula in the Nankai subduction zone. Potential physical mechanisms responsible for SSFs will have to explain this enigmatic dependence between duration and direction (Rubin, 2011). The ratio of along-strike SSFs propagating in the SSE direction to along-strike SSFs propagating in the opposite direction also seems to depend on the duration. While reversal SSFs are predominant in the 8 h time window, almost all SSFs propagate forward in the 32 h time window.
<figure><img src="image_39.png" /><figcaption>Fig. 2. Examples of SSFs detected in the LFE catalog for the different time windows of our analysis. Grey arrows show SSF direction determined by the algorithm.</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>Fig. 3. Examples of SSFs detected in the tremor catalog for the different time windows of our analysis. Grey arrows show SSF direction determined by the algorithm.</figcaption></figure>

The distribution of propagation velocity $V_{prop}$ is shown at the bottom of each subfigure in Figs. 4 and 5. In individual time windows, along-dip SSFs do not appear faster than along-strike SSFs. However, as the length $L$ of SSFs is limited by the extent of the fault, the average propagation velocity decreases with the duration $T$ (because $V_{prop} = L/T$). Accordingly, the limited number of distinct time windows in the analysis of the tremor catalog results in little difference in the propagation velocity (Fig. S1.f).

## 4. Physical characteristics of secondary slip fronts

Our detection algorithm provides a number of direct observations characterizing each SSF including start time, location, direction (Figs. 4, 5), duration $T$, length $L$, width $W$, pulse (i.e. actively slipping) length $l$ (Fig. 6), propagation velocity $V_{prop}$ (defined as $L/T$), and number of LFE or tremor detections $n_e$. The distributions of these SSF properties are shown in Fig. S1 (catalogs are available in the supplementary materials). Assessing other physical characteristics of SSFs, such as moment or stress drop, is challenging because these quantities require an estimate of the total slip associated with each SSF.

Individual SSFs are too small to generate a surface deformation signal detectable with GPS and hence most studies that have estimated slip for individual SSFs have relied on the total SSE moment and the fraction of LFE detections or bursts that occurred as part of the SSF (Rubin and Armbruster, 2013; Royer et al., 2015). Here, we adopt a similar approach and assume that the total moment $M_e$ associated with each LFE or tremor detection is the same (although different for LFEs and tremor). This assumption is demonstrably false for individual LFE repeats, which have been shown to have very different moments (Bostock et al.,

2015), and for individual tremor events, as their definition (i.e. the 2.5 min timescale used to separate individual events) is arbitrary and includes no contribution from waveform amplitudes (Wech and Creager, 2008). However, since moment estimates for each SSF include a minimum of 20 LFEs or tremors, it is likely that the average moment will be a more accurate measure of total SSF moment. Nevertheless, our estimates of the physical characteristic of SSFs (such as moment, slip, stress drop, slip rate) are tentative and caution should be exercised when interpreting them.

We obtain $M_e$, the average moment associated with each LFE or tremor, by dividing the moment $M_0$ (assessed by geodetic studies Schmidt and Gao, 2010; Krogstad and Schmidt, 2015) associated with an SSE by the number of events $N_e$ detected during the event:

$$M_e = \frac{M_0}{N_e} \qquad (1)$$

An implicit assumption behind equation (1) is that moment associated with aseismic slip that does not generate LFE or tremor detections is negligible. The moment, $M_{SSF}$, associated with one SSF, is defined as the sum of the moment associated with all the LFE or tremor detections that are part of that SSF,

$$M_{SSF} = \sum_{i=1}^{n_e} M_e(i), \qquad (2)$$

where $n_e$ is the number of LFE/tremor detections part of the given SSF. As we assumed that for all $i$, $M_e(i) = \frac{M_0}{N_e}$,

$$M_{SSF} = \frac{n_e}{N_e} M_0. \qquad (3)$$
<figure><img src="image_41.png" /><figcaption>Fig. 4. Detected secondary slip fronts in Cascadia from August 2009 to December 2015 for time windows of 4 h, 8 h, 16 h and 32 h (tremor catalog). White to gray color scale indicates the density of tremor detections. Arrows show lengths and directions of the SSFs propagating along strike (red) and along dip (black). Rose diagrams show the number of SSFs in the different directions. Histograms show velocity distributions for SSFs propagating along strike (red) and along dip (black). Red and black vertical lines superimposed on histograms show average velocity for SSFs propagating along strike and along dip, respectively.</figcaption></figure>

Note that possible deviations from constant $M_e(i)$ would tend to be averaged out over the number of events in equation (3) because $M_e(i)$ is sometimes larger and sometimes smaller than the assumed average $M_e = M_0/N_e$. We further estimate an average slip $d$ for each SSF,

$$d = \frac{M_{SSF}}{\mu A}, \qquad (4)$$

where $\mu$ is the shear modulus and $A = WL$ is the SSF slip area, with the SSF width $W$ and length $L$ both estimated from our analysis (Fig. 6). To be consistent with geodetic inversions, we use $\mu = 40$ GPa (Schmidt and Gao, 2010; Krogstad and Schmidt, 2015). The stress drop $\Delta\sigma$ is then given by

$$\Delta\sigma = \frac{4(\lambda + \mu)}{\pi(\lambda + 2\mu)} \mu \frac{d}{W}, \qquad (5)$$

where $\lambda$ is Lamé's parameter (Lay and Wallace, 1995), which we set to $\lambda = 40$ GPa. The total slip $d$ accumulates in the time the front propagates the distance $l$. Thus, we can deduce the slip rate $V_{slip}$ from the propagation velocity $V_{prop}$ as

$$V_{slip} = \frac{d}{l} V_{prop}. \qquad (6)$$

The distributions of assessed physical characteristics (moment, slip, stress drop and slip rate) are shown in Figs. S2 and S3 (catalogs in the supplementary materials). Because we did not estimate these parameters independently, their distributions are correlated with each other (Fig. S4). However, most of the correlation coefficients between the different parameters are relatively weak (i.e. below 0.6 in magnitude, see Fig. S4). As expected, high correlation coefficients are found for slip and stress drop, since the stress drop is estimated from slip (eq. (5)), and between the different spatial dimensions ($L, W$ and $l$), implying that the aspect ratios of SSFs do not vary appreciably. Significant anticorrelations are observed between propagation rate $V_{prop}$ and duration $T$, as well as between the spatial dimensions and both slip $d$ and slip rate $V_{slip}$. It is also noteworthy that stress drop $\Delta\sigma$ appears to be independent of moment $M$ in both catalogs (Fig. S4).

The median stress drop among all SSFs is 5.8 kPa, which is very consistent with stress drop estimates obtained for RTRs using strainmeters (8 kPa) (Hawthorne et al., 2016). These values of stress drop are several times smaller than typical stress drops in SSEs (10–100 kPa), and over 2 orders of magnitude smaller than regular earthquakes (Gao et al., 2012). The estimated stress drops are also all consistently smaller than effective normal stress (0.2–3 MPa) inferred in LFE source regions (Liu and Rice, 2007; Hawthorne and Rubin, 2013; Royer et al., 2015). The median slip
<figure><img src="image_42.png" /><figcaption>**Fig. 5.** Detected secondary slip fronts during the 2003, 2004 and 2005 SSEs for time windows of 30 min, 1 h, 2, h, 4 h, 8 h, 16 h and 32 h (LFE catalog). White to gray color scale indicates the relative density of LFEs. Blue thick arrow shows the approximate direction of the three SSEs (45° west from north). Red arrows show lengths and directions of SSFs propagating in the direction of SSE propagation ±45° (along strike). Black arrows show lengths and directions of SSFs propagating perpendicular to the direction of SSE propagation ±45° (along dip). Rose diagrams show the number of SSFs in the different directions. Histograms show velocity distributions for SSFs propagating along strike (red) and along dip (black). Red and black vertical lines superimposed on histograms show average velocity for SSFs propagating along strike and along dip, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>**Fig. 6.** Scheme illustrating how SSF length (*L*), width (*W*), and pulse length (*l*) are deduced from the analysis. *L* is defined as the length of the best projection axis (parallel to the propagation direction) (a and b). *l* is defined as the 2σ dispersion from the linear regression of the distance along the best projection axis as a function of time (a). *W* is defined as the 2σ spatial dispersion (perpendicular to the propagation direction) from the best projection axis (b).</figcaption></figure>

is 3 mm and the median slip rate is 1.1 mm/h, which is much faster than slip rates inferred for SSEs (on the order of 1–2 mm/day Wech and Bartlow, 2014). The observation that SSF slip rates are significantly faster than SSE slip rates is consistent with previous studies (Rubin and Armbruster, 2013; Royer et al., 2015; Hawthorne et al., 2016) and implies that the fault rheology limits the main front speed and then allows secondary fronts to slip faster. As acknowledged in the introduction, such a time dependent behavior is not easy to reproduce with a single state variable

friction law (Rubin, 2011). SSF slip rates may then be used to constrain physical models of slow slip phenomena.

Though variations in the moments and slip distributions of SSFs range within a single order of magnitude (4.10¹⁶ to 4.10¹⁷ N.m and 0.1 to 40 mm) for all SSFs, stress drop varies by over 3 orders of magnitude (from 0.1 to 300 kPa), as does slip rate (between 0.02 and 30 mm/h) (Figs. S2 and S3). We recall here that our estimates of stress drop are based on the original assumption that SSF moment scales with the number of LFE/tremor detections, which could be inaccurate for some SSFs. Large scatter in slip, stress drop and slip rate is also expected from equations (4), (5) and (6) as a result of multiplying quantities ($M_{SSF}$, *L*, *W*, *l* and $V_{prop}$) that have large variability. Therefore a large part of the scatter in slip rate and stress drop might be attributed to assumption inaccuracy and propagation of uncertainties, but some of this scatter could also be sourced to changes in fault properties. As SSF areas widely overlap — this is particularly clear in the LFE catalog (Fig. 5) — the varying SSF properties might suggest that, in addition to being able to slip at different rates, the same fault area is also able to release very different amounts of stress.

## 5. Moment-duration scaling

Ide et al. (2007) suggested that slow slip phenomena — including LFEs, very low frequency earthquakes (VLFEs) and SSEs — follow a moment-duration scaling law distinct from regu-
<figure><img src="image_44.png" /><figcaption>**Fig. 7.** Moment-duration scaling law for regular earthquakes and slow slip phenomena (gray shaded areas). Colored solid lines are isolines of equal moment rate. Black stars indicate regular earthquakes over a large range of magnitude (compilation of Peng and Gomberg, 2010 augmented with M4.8–5.8 events from the GCMT catalog and recent mega-earthquakes, Delouis et al., 2010; Bletery et al., 2014; Bletery et al., 2016). Gray stars indicate tsunami earthquakes (compilation from Peng and Gomberg, 2010). Blue dots are SSEs (compilation of Peng and Gomberg, 2010 augmented with events from Gao et al., 2012 and Liu, 2014). Orange dots are VLFEs (Ito and Obara, 2006; Matsuzawa et al., 2009; Ide and Yabe, 2014). Orange circles are unnamed 20–200 s long duration earthquakes assimilated to VLFEs (Ide et al., 2008). Red circles are LFE moment and duration averaged over 100 (Thomas et al., 2016) and >100,000 events (Bostock et al., 2015). Cyan dashed line shows an approximate moment threshold ($M_w \approx 5.7–6.0$) for geodetically detectable events. Red dashed line shows the limit of seismically detectable events (see text). In agreement with moment of large SSFs determined from strain meters (green square, Hawthorne et al., 2016), SSFs detected in this study (green dots) span a moment-duration region that cannot be directly recorded either by common seismic or geodetic instrumentation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

lar earthquakes. While regular (fast) earthquakes have moment proportional to the cube of their duration, slow slip phenomena would have moment proportional to duration (with $M/T = 10^{12}–10^{13}$ N.m/s) (Fig. 7). In this representation, important gaps separate the various slow slip phenomena (Ide et al., 2007). The question of whether these gaps result from detection limitations or from the inherent behavior of fault slip is of importance because it could provide us with physical constraints on earthquake processes (Ide et al., 2007; Ide, 2008; Ide et al., 2008; Peng and Gomberg, 2010; Gao et al., 2012; Ide and Yabe, 2014; Liu, 2014; Gomberg et al., 2016).

In the far field, seismic amplitudes scale with moment rate (Aki and Richards, 2002). Thus, the threshold for seismically detectable signals follows the noise level in the different frequency bands. We use the USGS low noise model (Peterson et al., 1993; Ide, 2014) and empirically calibrate it on the smallest moment rate events that have been seismically detected (VLFEs from Ide and Yabe, 2014) to determine a threshold for seismically detectable events (red dashed line in Fig. 7). Existing geodetic instrumentation has proven capable of detecting magnitude 5.7–6.0 events and above regardless of their duration (Peng and Gomberg, 2010) (blue dashed line in Fig. 7). The gap between LFEs and VLFEs (Fig. 7) may be attributed to particularly large noise in the frequency band corresponding to the gap (Ide et al., 2008) while the gap between the VLFEs and SSEs corresponds to moments too low to be detected by geodesy (Peng and Gomberg, 2010) and characteristic periods too long to be detected in seismograms (Ide et al., 2008).

The large range of duration and moment among the SSFs we detected follows the slow slip scaling law and fills the gap between

VLFEs and SSEs, where signal is not directly detectable either from seismic or geodetic records (Fig. 7). The vertical (moment) threshold in SSFs corresponds to the minimum number of LFE/tremor detections (20) in one SSF times $M_e$ and the distinct horizontal lines in SSF duration correspond to our different time windows. These two noncontinuous characteristics are consequences of our approach. Interestingly, SSFs appear to follow the same threshold of constant moment rate — on the order of $10^{13}$ N.m/s — as the other slow slip phenomena (Fig. 7). The fact that SSFs follow the slow slip phenomena scaling law may be a consequence of estimating SSF moment from moment of SSEs (eq. (3)), as SSEs follow this law. Nevertheless, it is worth noting that our moment estimates are in agreement with SSF moments estimated from strain meters (Hawthorne et al., 2016) (green square in Fig. 7).

While individual slow slip phenomena do not appear to have moment proportional to duration (e.g. Ide et al., 2008; Bostock et al., 2015; Gomberg et al., 2016), they all fit a constant moment rate scaling law (within less than 2 orders of magnitude variation) when considered together on a 10 orders of magnitude scale (Fig. 7). More generally, fast and slow slip phenomena appear to be well distinguished in Fig. 7. Fault slip processes that would range between the two scaling laws are, in theory, detectable — and should therefore have been detected — by seismology, geodesy or both (Fig. 7). The large gap in detection of fault slip processes between the two scaling laws suggests that earthquakes and slow slip phenomena are two very distinct fault slip processes.

## 6. Conclusion

We identified a total of 1076 SSFs in Cascadia at timescales between 30 min and 32 h. We cataloged their duration, width, length, pulse length and time of occurrence and provided estimates for some of their physical properties such as moment, slip, slip rate, and stress drop. The median moment is $5.3 \times 10^{16}$ N.m (corresponding to a moment magnitude of 5.1), the median slip is 3 mm, the median slip rate is 1.1 mm/h and the median stress drop is 5.8 kPa. Moment estimates suggest that SSFs bridge a major gap between seismically and geodetically detectable events in the slow slip phenomena scaling law (Ide et al., 2007). As previously observed with different approaches (Ghosh et al., 2010; Houston et al., 2011; Obara et al., 2012), we found that short duration SSFs dominantly propagate along dip while long duration SSFs dominantly propagate along strike. We also found that inferred SSF stress drops and slip rates range over 3 orders of magnitude, suggesting that the same fault region could be able to slip at very different rates and to release very different stresses. These findings may be used to constrain models of slow earthquake processes.

## Acknowledgements

This work was funded by the NSF grant EAR-1520238. The tremor catalog was made available online by A. Wech and K. Crea-ger (tunk.ess.washington.edu/map_display/). We thank G. Savard and M. Bostock for providing the LFE catalog.

## Appendix A. Supplementary material

Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.epsl.2017.01.046.

## References

* Aki, K., Richards, P.G., 2002. Quantitative Seismology. Vol. 1.

* Bartlow, N.M., Miyazaki, S., Bradley, A.M., Segall, P., 2011. Space-time correlation of slip and tremor during the 2009 Cascadia slow slip event. Geophys. Res. Lett. 38 (18).
Bletery, Q., Sladen, A., Delois, B., Vallée, M., Nocquet, J.-M., Rolland, L., Jiang, J., 2014. A detailed source model for the Mw9.0 Tohoku-Oki earthquake reconciling geodesy, seismology, and tsunami records. *J. Geophys. Res., Solid Earth* 119 (10), 7636–7653.

Bletery, Q., Sladen, A., Jiang, J., Simons, M., 2016. A Bayesian source model for the 2004 great Sumatra-Andaman earthquake. *J. Geophys. Res., Solid Earth* 121 (7), 5116–5135.

Bostock, M., Thomas, A., Savard, G., Chuang, L., Rubin, A., 2015. Magnitudes and moment-duration scaling of low-frequency earthquakes beneath southern Vancouver island. *J. Geophys. Res., Solid Earth* 120 (9), 6329–6350.

Chiu, S.L., 1994. Fuzzy model identification based on cluster estimation. *J. Intell. Fuzzy Syst.* 2 (3), 267–278.

Delois, B., Nocquet, J.-M., Vallée, M., 2010. Slip distribution of the February 27, 2010 Mw = 8.8 Maule Earthquake, central Chile, from static and high-rate GPS, InSAR, and broadband teleseismic data. *Geophys. Res. Lett.* 37 (17).

Dragert, H., Wang, K., James, T.S., 2001. A silent slip event on the deeper Cascadia subduction interface. *Science* 292 (5521), 1525–1528.

Frank, W.B., Shapiro, N.M., Husker, A.L., Kostogllovov, V., Romanenko, A., Campillo, M., 2014. Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico. *J. Geophys. Res., Solid Earth* 119 (10), 7686–7700.

Gao, H., Schmidt, D.A., Weldon, R.J., 2012. Scaling relationships of source parameters for slow slip events. *Bull. Seismol. Soc. Am.* 102 (1), 352–360.

Ghosh, A., Vidale, J.E., Sweet, J.R., Creager, K.C., Wech, A.G., Houston, H., Brodsky, E.E., 2010. Rapid, continuous streaking of tremor in Cascadia. *Geochem. Geophys. Geosyst.* 11 (12).

Gomberg, J., Wech, A., Creager, K., Obara, K., Agnew, D., 2016. Reconsidering earthquake scaling. *Geophys. Res. Lett.* 43 (12), 6243–6251.

Hawthorne, J.C., Bostock, M.G., Royer, A.A., Thomas, A.M., 2016. Variations in slow slip moment rate associated with rapid tremor reversals in Cascadia. *Geochem. Geophys. Geosyst.* 17 (12), 4899–4919.

Hawthorne, J., Rubin, A., 2013. Tidal modulation and back-propagating fronts in slow slip events simulated with a velocity-weakening to velocity-strengthening friction law. *J. Geophys. Res., Solid Earth* 118 (3), 1216–1239.

Houston, H., 2015. Low friction and fault weakening revealed by rising sensitivity of tremor to tidal stress. *Nat. Geosci.* 8 (5), 409–415.

Houston, H., Delbridge, B.G., Wech, A.G., Creager, K.C., 2011. Rapid tremor reversals in Cascadia generated by a weakened plate interface. *Nat. Geosci.* 4 (6), 404–409.

Ide, S., 2008. A Brownian walk model for slow earthquakes. *Geophys. Res. Lett.* 35 (17).

Ide, S., 2010. Striations, duration, migration and tidal response in deep tremor. Nature 466 (7304), 356–359.

Ide, S., 2014. Modeling fast and slow earthquakes at various scales. *Proc. Jpn. Acad., Ser. B* 90 (8), 259–277.

Ide, S., Beroza, G.C., Shelly, D.R., Uchide, T., 2007. A scaling law for slow earthquakes. *Nature* 447 (7140), 76–79.

Ide, S., Imanishi, K., Yoshida, Y., Beroza, G.C., Shelly, D.R., 2008. Bridging the gap between seismically and geodetically detected slow earthquakes. *Geophys. Res. Lett.* 35 (10).

Ide, S., Yabe, S., 2014. Universality of slow earthquakes in the very low frequency band. *Geophys. Res. Lett.* 41 (8), 2786–2793.

Ito, Y., Obara, K., 2006. Very low frequency earthquakes within accretionary prisms are very low stress-drop earthquakes. *Geophys. Res. Lett.* 33 (9).

Kao, H., Shan, S.-J., Dragert, H., Rogers, G., Cassidy, J.F., Wang, K., James, T.S., Ramachandran, K., 2006. Spatial-temporal patterns of seismic tremors in northern Cascadia. *J. Geophys. Res., Solid Earth* 111 (B3).

Krogstad, R.D., Schmidt, D.A., 2015. Assessing the updpi spatial offset of tremor and slip during ETS events in Cascadia. In: 2015 AGU Fall Meeting. AGU.

Lay, T., Wallace, T.C., 1995. Modern Global Seismology. International Geophysics Series, vol. 58. Academic Press.

Liu, Y., 2014. Source scaling relations and along-strike segmentation of slow slip events in a 3-d subduction fault model. *J. Geophys. Res., Solid Earth* 119 (8), 6512–6533.

Liu, Y., Rice, J.R., 2007. Spontaneous and triggered aseismic deformation transients in a subduction fault model. *J. Geophys. Res., Solid Earth* 112 (B9).

Matsuzawa, T., Obara, K., Maeda, T., 2009. Source duration of deep very low frequency earthquakes in western Shikoku, Japan. *J. Geophys. Res., Solid Earth* 114 (B11).

Obara, K., 2002. Nonvolcanic deep tremor associated with subduction in southwest Japan. *Science* 296 (5573), 1679–1681.

Obara, K., Matsuzawa, T., Tanaka, S., Maeda, T., 2012. Depth-dependent mode of tremor migration beneath Kii peninsula, Nankai subduction zone. *Geophys. Res. Lett.* 39 (10).

Peng, Z., Gomberg, J., 2010. An integrated perspective of the continuum between earthquakes and slow-slip phenomena. *Nat. Geosci.* 3 (9), 599–607.

Peng, Y., Rubin, A.M., 2016. High-resolution images of tremor migrations beneath the Olympic Peninsula from stacked array of arrays seismic data. *Geochem. Geophys. Geosyst.* 17, 587–601.

Peng, Y., Rubin, A.M., 2017. Intermittent tremor migrations beneath Guerrero, Mexico and implications for fault healing within the slow slip zone. *Geophys. Res. Lett.*

Peng, Y., Rubin, A.M., Bostock, M.G., Armbruster, J.G., 2015a. High-resolution imaging of rapid tremor migrations beneath southern Vancouver island using cross-station cross correlations. *J. Geophys. Res., Solid Earth* 120 (6), 4317–4332.

Peng, Z., Shelly, D.R., Ellsworth, W.L., 2015b. Delayed dynamic triggering of deep tremor along the Parkfield-Cholame section of the San Andreas Fault following the 2014 M6.0 South Napa earthquake. *Geophys. Res. Lett.* 42 (19), 7916–7922.

Peterson, J., et al., 1993. Observations and Modeling of Seismic Background Noise. *US Geological Survey, Albuquerque, New Mexico*.

Radiguet, M., Cotton, F., Vergnolle, M., Campillo, M., Walpersdorf, A., Cotte, N., Kostogllovov, V., 2012. Slow slip events and strain accumulation in the Guerrero gap, Mexico. *J. Geophys. Res., Solid Earth* 117 (B4).

Royer, A., Thomas, A., Bostock, M., 2015. Tidal modulation and triggering of low-frequency earthquakes in northern Cascadia. *J. Geophys. Res., Solid Earth* 120 (1), 384–405.

Rubin, A.M., 2011. Designer friction laws for bimodal slow slip propagation speeds. *Geochem. Geophys. Geosyst.* 12 (4).

Rubin, A.M., Armbruster, J.G., 2013. Imaging slow slip fronts in Cascadia with high precision cross-station tremor locations. *Geochem. Geophys. Geosyst.* 14 (12), 5371–5392.

Ryberg, T., Haberland, C., Fuis, G., Ellsworth, W., Shelly, D., 2010. Locating non-volcanic tremor along the San Andreas Fault using a multiple array source imaging technique. *Geophys. J. Int.* 183 (3), 1485–1500.

Savard, G., Bostock, M.G., 2015. Detection and location of low-frequency earthquakes using cross-station correlation. *Bull. Seismol. Soc. Am.*.

Schmidt, D., Gao, H., 2010. Source parameters and time-dependent slip distributions of slow slip events on the Cascadia subduction zone from 1998 to 2008. *J. Geophys. Res., Solid Earth* 115 (B4).

Shelly, D.R., 2009. Possible deep fault slip preceding the 2004 Parkfield earthquake, inferred from detailed observations of tectonic tremor. *Geophys. Res. Lett.* 36 (17).

Shelly, D.R., 2010. Migrating tremors illuminate complex deformation beneath the seismogenic San Andreas Fault. *Nature* 463 (7281), 648–652.

Shelly, D.R., Beroza, G.C., Ide, S., 2007a. Complex evolution of transient slip derived from precise tremor locations in western Shikoku, Japan. *Geochem. Geophys. Geosyst.* 8 (10).

Shelly, D.R., Beroza, G.C., Ide, S., 2007b. Non-volcanic tremor and low-frequency earthquake swarms. *Nature* 446 (7133), 305–307.

Thomas, A.M., Beroza, G.C., Shelly, D.R., 2016. Constraints on the source parameters of low-frequency earthquakes on the San Andreas Fault. *Geophys. Res. Lett.* 41 (2), 392–398.

Thomas, T.W., Vidale, J.E., Houston, H., Creager, K.C., Sweet, J.R., Ghosh, A., 2013. Evidence for tidal triggering of high-amplitude rapid tremor reversals and tremor streaks in northern Cascadia. *Geophys. Res. Lett.* 40 (16), 4254–4259.

Ueno, T., Maeda, T., Obara, K., Asano, Y., Takeda, T., 2010. Migration of low-frequency tremors revealed from multiple-array analyses in western Shikoku, Japan. *J. Geophys. Res., Solid Earth* 115 (B9).

Wech, A.G., 2010. Interactive tremor monitoring. *Seismol. Res. Lett.* 81 (4), 664–669.

Wech, A.G., Bartlow, N.M., 2014. Slip rate and tremor genesis in Cascadia. *Geophys. Res. Lett.* 41 (2), 392–398.

Wech, A.G., Creager, K.C., 2008. Automated detection and location of Cascadia tremor. *Geophys. Res. Lett.* 35 (20).

Yager, R.R., Filev, D.P., 1994. Generation of fuzzy rules by mountain clustering. *J. Intell. Fuzzy Syst.* 2 (3), 209–219.
## 3 La coalescence de séismes lents comme possible mécanisme d'initiation des grands séismes

Lors de mes travaux sur les tremors à l'Université d'Oregon, j'ai été frappé par la distribution des tremors associés à un séisme lent particulier, survenu en 2013 au Nord de la faille des Cascades (évènement entouré en noir sur la Figure III.1). Alors que les séismes lents se propagent normalement, soit unilatéralement, soit bilatéralement, formant dans la distribution des tremors représentée dans un diagramme latitude-temps, soit une droite, soit un éventail ouvert vers la droite, les tremors associés à cet évènement forment un éventail vers la gauche, comme si ce dernier défiait la causalité (Figure III.1).

<figure><img src="image_45.png" /><figcaption>**Figure III.1** – Latitude des tremors en fonction du temps le long de la faille des Cascades [Wech and Creager, 2008; Wech, 2010]. L'ellipse noire indique l'évènement particulier de 2013.</figcaption></figure>

Après mon recrutement à Géoazur, j'ai voulu savoir si le glissement lent associé à cet épi-sodes de tremors présentait également cette étrange distribution. J'ai, pour cela, utilisé les outils développés par Jean-Mathieu Nocquet pour inverser le glissement journalier le long de la faille. Nous avons trouvé, dans l'étude qui suit, que le glissement lent présentait effectivement une distribution exceptionnelle, originant non pas en 2 mais en 3 sous évènements distincts coalesçant
### III.3 La coalescence de séismes lents comme possible mécanisme d'initiation des grands séismes

en 1 seul. De manière encore plus intéressante, nous avons observé que lors des 2 coalescences successives de ces sous-évènements (observées de manière concordante dans l'inversion GNSS et dans la distribution des tremors), l'énergie (le moment) relâchée a significativement augmenté, suggérant que la coalescence de séismes lents puisse être un mécanisme d'initiation des grands séismes.

Cette étude a été réalisée à Géoazur, en collaboration avec Jean-Mathieu Nocquet. Elle a été publiée en 2020 dans le journal Nature Communications [Bletery and Nocquet, 2020].
<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

Check for updates

ARTICLE

https://doi.org/10.1038/s41467-020-15494-4

OPEN

# Slip bursts during coalescence of slow slip events in Cascadia

Quentin Bletery<fnref n="1" />✉ & Jean-Mathieu Nocquet<fnref n="1" />,²

Both laboratory experiments and dynamic simulations suggest that earthquakes can be preceded by a precursory phase of slow slip. Observing processes leading to an acceleration or spreading of slow slip along faults is therefore key to understand the dynamics potentially leading to seismic ruptures. Here, we use continuous GPS measurements of the ground displacement to image the daily slip along the fault beneath Vancouver Island during a slow slip event in 2013. We image the coalescence of three originally distinct slow slip fronts merging together. We show that during coalescence phases lasting for 2 to 5 days, the rate of energy (moment) release significantly increases. This observation supports the view proposed by theoretical and experimental studies that the coalescence of slow slip fronts is a possible mechanism for initiating earthquakes.

<footnote n="1">¹Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Géoazur, 250 rue Albert Einstein, 06560 Valbonne, France. ²Institut de Physique du Globe de Paris, Université de Paris, CNRS, 75238 Paris, France. ✉email: bletery@geoazur.unice.fr</footnote>
In western North America, the northern Cascadia subduction fault experiences slow slip events (SSEs) approximately every 14 months¹. As slow slip develops, it generates a particular type of micro-earthquakes, known as low-frequency earthquakes that superpose as tectonic tremors¹⁻³. The space-time correlation between slow slip and tremors⁴,⁵ gives the opportunity to study the details of the evolution of SSEs, taking advantage of independent measurements of the ground displacement and tremor activity⁶⁻¹². Here we focus on a complex SSE that occurred beneath Vancouver Island in September–October 2013. We use Global Positioning System (GPS) time series of the ground displacement to image the daily evolution of slip along the Cascadia subduction fault (Figs. 1 and 2) in relation to the daily distribution of tremors cataloged by the Pacific Northwest Seismic Network (PNSN)¹³,¹⁴.

## Results

**Two coalescence episodes.** We find that slow slip initiates in the area of Seattle on September 7, 2013, associated with intense tremor activity (Supplementary Fig. 3). During the first 2 weeks, most of the slip and tremors remain confined to the same 60 × 30 km² area, with a maximum activity on September 22 (Supplementary Figs. 3 and 4). Starting on September 15, tremors also emerge ~250 km from that area at latitude 50° under Vancouver Island with low-amplitude slip sparsely detected from September 19 onward (Supplementary Fig. 4). Tremor locations indicate a southeastward, along-strike, and constant depth migration at ~10 km per day. On September 26, slip catches up with tremors and extends southeastward. On September 24, an additional area of tremors arises between the two SSE areas, in the southern part of Vancouver island, soon followed by slow slip. On September

<figure><img src="image_48.png" /><figcaption>**Fig. 1** First coalescence episode. Color maps show the slow slip rate inverted along the fault on each day from September 25 to September 28, 2013. Blue dots show the location of tremors recorded on the corresponding day. Isolines indicate the fault depth from 20 to 80 km.</figcaption></figure>
<figure><img src="image_49.png" /><figcaption>**Fig. 2** Second coalescence episode. Color maps show the slow slip rate inverted along the fault on each day from September 29 to October 4, 2013. Blue dots show the location of tremors recorded on the corresponding day. Isolines indicate the fault depth from 20 to 80 km.</figcaption></figure>

26, we see three clearly distinct patches of slip, each of them associated with a tremor swarm (Fig. 1). From September 27, the region located in between the southern and central patches starts slipping. On September 28, the two southern patches are indistinguishable one from the other (Fig. 1). During this merging phase, the moment rate release raises by a factor 10 in 3 days and the daily maximum slip increases by a factor 3 in 2 days (Fig. 3a). As an independent check for this acceleration, we find that the average daily incremental displacement recorded by the GPS stations close to the merging area is multiplied by a factor 8 (Fig. 3b) (see “Methods—Edge filter”). An interesting aspect is that, before the merging phase, the slow slip at the southern patch was fading out (see decrease in moment rate release after September 22, Fig. 3a), suggesting that the coalescence of the two previously distinct slipping areas is the actual cause for the observed increase in moment rate release. Following the first merging, the moment release slightly diminishes (Figs. 2 and 3a). A second merging phase occurs quickly afterwards (Fig. 2). This latter merging phase is also associated with a significant pulse in the moment release rate (2× between October 1 and October 4),

which eventually decreases to end the sequence on October 12 (Fig. 3, Supplementary Fig. 5).

**Interaction between distant slip fronts?** The 2013 Fall sequence highlights two successive coalescence episodes of pre-existing distinct slipping areas. Though jumping and halting are common in SSEs⁶, merging episodes similar to this sequence have not been observed during other SSEs in Cascadia or elsewhere. Nevertheless, the coalescence of nearby slipping areas might not be uncommon. In a rate-and-state framework¹⁵,¹⁶, as slip propagates from one area toward another, it progressively meets increasing positive shear stress generated by the counterpart patch, which in turn enhances its slip and promotes further propagation, until both slipping areas coalesce. For both merging phases during the 2013 sequence, our inversion suggests that coalescence initiates when the distance between the two SSE areas is on the order of one length-scale of their size. However, interaction at larger distance cannot be ruled out since tremors —and by extension SSEs—are known to be sensitive to extremely small stress
<figure><img src="image_50.png" /></figure>

**Fig. 3** Increase in moment release during both coalescence episodes. **a** Evolution of the moment rate release during the entire sequence. Gray histograms indicate the relative number of tremors on each day. The red thin curve shows the maximum daily slip on the fault every day. The blue thick curve shows the moment rate release. **b** Average daily displacement recorded on GPS stations close to (blue curve) and away from (red curve) the merging area after application of an edge filter, highlighting an increase in GPS data amplitude within the stations located close to the coalescence episodes. Stations used in both stacks are shown in Supplementary Fig. 1 (blue and red stations close to and away from the merging area, respectively).

perturbations, such as those generated by tides¹⁷. The southern migration of the northern slip front or the onset of the center one possibly reflects such long-distance interaction.

**Enhanced slip during coalescence episodes.** Both merging phases are associated with a burst in moment rate release resulting from a combination of a growth of the slipping area filling the space between the two patches and an increase in daily slip velocity (Fig. 3a). Moment release bursts result in an acceleration of the daily incremental displacement at GPS sites located above the merging area, directly seen in the GPS data (Fig. 3b). The period of merging between September 27 and October 5 (9 days) contributes to more than half (59%) of the total moment released during the entire sequence (47 days). Furthermore, peak moment rates on September 28 and October 4 represent a third (33%) of the moment released during the merging period. These ratios suggest a merging-induced burst-like behavior of slip at a daily time scale. Slip might, however, experience sub-daily variations. Indeed, analyses of seismic data indicate that microseismic activity, presumably driven by slow slip, is clustered in minute-to-hour-long bursts¹⁸,¹⁹. The daily incremental slip inferred here is likely an average of multiple shorter episodic events, with the implication that slip rates may punctually be significantly faster than the daily slip velocities inferred from our inversion.

Discussion

Laboratory-derived constitutive friction laws predict that, for certain combinations of parameters referred to as velocity-weakening regimes, when sliding initiates, the friction drops proportionally to the logarithm of the sliding velocity, which in turn leads to enhanced sliding velocity¹⁵,¹⁶. Applied to the earthquake initiation problem, dynamic simulations have shown that the conditions required for slip to become unstable are functions—among other parameters—of the slip rate or stress perturbation in the nucleation area¹⁵,¹⁶. The coalescence of slow slip fronts therefore provides a mechanism counteracting the usual damping of slip observed during SSEs. This mechanism promotes instability on the fault through a slip acceleration and a stress increase at the tip of the locked area, which is larger than the summed contribution of individual fronts.

Based on laboratory experiments and numerical simulation, Dieterich intuited 40 years ago that slip velocity must increase if two previously independent zones of slip approach and coalesce to form a single zone²⁰. His reasoning was that, because stress on a dislocation is controlled by the ratio of slip over the length of the slipping area, the rate of slip must accelerate to double the displacement as the segments merge to form a single zone. Dieterich²⁰ and more recent studies²¹,²² proposed the coalescence of slow slip fronts as a possible mechanism to initiate earthquakes. While not the case here, the mechanism might have been at play during SSEs that have been proposed to trigger earthquakes.
Among the examples of recent large earthquakes that have been preceded by SSEs<sup>23-25</sup>, multiple slow slips accompanied by intense seismicity have been documented from months to days before the M<sub>w</sub> 8.1 April 1, 2014 Iquique earthquake<sup>26,27</sup>. Although the merging of SSEs has not been geodetically identified, the sequence shows some similarities with the 2013 Cascadia sequence. The migration of small earthquakes and activity of repeating earthquakes both suggest the existence of two slow slip areas separated by ~120 km in January and February 2014 north and south of the epicenter. While the northern area appears to be stationary and shows repeating earthquake activity until March 23, clustered and intense seismicity following the March 16 M<sub>w</sub> 6.7 earthquake migrates from south to north<sup>27</sup>, leaving the possibility for the two slow slip areas to have merged near the future epicenter.

The 2013 Cascadia sequence provides an observational evidence of a phenomenon previously proposed by physical models<sup>20-22</sup>. It demonstrates that periods of SSE merging are key periods during which slip is enhanced, possibly bringing the fault closer to seismic failure. In the quest for earthquake precursors, it provides a clue on where to pay attention.

## Methods

**GPS data.** We use 59 GPS daily time series from the final Geodesy Advancing Geosciences and EarthScope combined solution<sup>28</sup> available from the UNAVCO data center and added 3 stations (ELIZ, WOST, GLDR) from the Pacific Northwest Geodetic Array of the Central Washington University, north of Vancouver island where coverage is sparse. The location of the 62 used stations is shown in Supplementary Fig. 1. We remove outliers and seasonal signals over the whole time series. We detrend the time series, so that the velocity during the weeks before and after the 2013 event is zero. This procedure ensures that the transient displacement recorded at the GPS sites reflects positive reverse slip at the plate interface. Regional common mode motion evaluated by stacking time series for 10 sites with complete data and good repeatability located 100–300 km away for the slip area was found to be negligible. The obtained time series are shown in Supplementary Figs. 7–17 and discussed in Supplementary Note 1.

**Green's functions.** We use a curved geometry for the subduction interface based on the slab2.0 model for Cascadia<sup>29</sup> discretized into 307 quasi-equilateral triangular subfaults, ranging from 46°N to 50°N in latitude and from the trench down to 50 km. The rake is calculated on each subfault accounting for the rotation of the North American plate with respect to the Juan de Fuca plate (Euler rotation pole: −111.7°E, 32.0°N<sup>30</sup>). We calculate the transfer matrix g relating unit slip in the rake direction at each triangular subfault to the displacement components at the GPS sites using the solution for triangular dislocation elements in a uniform elastic half space<sup>31</sup>.

**Slip time-dependent inversion.** We invert for daily incremental slip from August 27, 2013 to October 12, 2013 using only the horizontal components. Our approach relates the vector of cumulative displacements $d_k$ at GPS sites since August 27 at date k to the sum of all daily incremental slip $\delta s_i$ before date k multiplied by the Green's function g:

$$ \sum_{i=1}^{k} g \delta s_i = d_k \quad (1) $$

Repeated for all dates k, Eq. (1) leads to a linear system $G\delta s = d$, where G is a block-triangular matrix made of g, d and $\delta s$ are the displacement vector and modeled daily slip resulting from concatenating all $d_k$ and $\delta s_i$, respectively. The chosen formulation allows to solve consistently for the daily slip at every subfault in a single inversion. Two important aspects of our inversion procedure are a non-negativity constraint<sup>32</sup> that prescribes backward slip increment at all time steps, strongly reducing the solution space<sup>33</sup> and the absence of temporal smoothing that allows a better restitution of slip acceleration. Spatial regularization constraints are imposed by the mean of a model covariance matrix controlling the level of damping and smoothing with respect to an a priori model<sup>8,34,35</sup> here taken as 0. The model covariance matrix is taken as an isotropic decreasing exponential<sup>8</sup>:

$$ C_{m,i,j} = \left( \sigma \frac{d_0}{D_c} \right)^2 \exp(-d_{i,j}/D_c) \quad (2) $$

where $C_{m,i,j}$ are the model covariance matrix elements corresponding to the i-th and j-th subfaults, $\sigma$ a constant (unit in mm day<sup>-1/2</sup>) controlling the weight of the regularization, $d_{i,j}$ the distance between the center of subfaults i and j, and $D_c$ is a correlation length controlling the level of spatial smoothing<sup>36</sup>. $d_0$ is a reference

distance taken as the mean length of the triangles. No temporal smoothing is added here aside from the non-negativity constraint imposing that slip must grow through time. We used $D_c = 50$ km, a value on the order of the distance between the center of two adjacent subfaults. We tested different values of $\sigma$ and plotted the gain in misfit reduction as a function of regularization ($\sigma$). We obtained a figure commonly known as an L-curve (Supplementary Fig. 2). We chose $\sigma = 3$ mm day<sup>-1/2</sup>, the point at the corner of the L-curve for which the misfit minimization becomes marginal when regularization increases.

**Tremor data.** We analyze the obtained slip distribution in light of the tremor activity (Figs. 1–3, Supplementary Figs. 3–6). We use the tremor catalog from the PNSN<sup>13,14</sup>. Because GPS data are daily averages of the positions measured from 00:00:00 to 23:59:30 UTC, they can be seen as the positions at 12:00:00 UTC every day. Therefore, the daily slip on the fault we image on day i should be seen as the slip from day (i − 1) at 12:00:00 to day i at 11:59:59 (Figs. 1 and 2, Supplementary Figs. 3–5). Accordingly, the tremors represented on the subfigure of day i are tremors recorded from day (i − 1) at 12:00:00 to day i at 11:59:59 (Figs. 1–3, Supplementary Figs. 3–5).

**Moment rate release.** The moment rate release (Fig. 3a) is calculated by summing the contributions of all subfaults assuming a uniform elastic Earth with an elastic shear modulus of 30 GPa. The cumulative moment is $1.10 \times 10^{19}$ N m, equivalent to a moment magnitude $M_w = 6.6$, the maximum cumulative slip is 35.3 mm, the maximum slip rate $2.7$ mm day<sup>-1</sup>. Because the dimensions of all the subfaults are equal, the moment rate release is strictly proportional to the average slip rate integrated over the ruptured surface. The daily number of tremors correlates fairly well with the moment release rate (coefficient of correlation: 0.63) and even better with daily maximum slip rate (coefficient of correlation: 0.76) (Fig. 3a). Furthermore, the two merging phases are associated with spikes in tremor activity (Figs. 1–3) and the day-to-day spatial correlation between slow slip and the independently inferred distribution of tremors (Figs. 1 and 2, Supplementary Figs. 3–5, Supplementary Movie 1) suggests that the main features of our inverted slip distribution are robust. Additional resolution tests can be found in Supplementary Information in the online version of the article (Supplementary Note 2, Supplementary Figs. 18–22).

**Edge filter.** In order to verify that the moment rate increase found in our inversion does not arise from fitting noise, we investigate the high-frequency properties of our time series. We apply a one-dimensional edge filter<sup>37</sup> to the north and east components of the GPS displacement time series. Edge filtering is a technique widely used in image processing that allows to retrieve sharp gradients in noisy data. We differentiate the filtered time series, take the norm, stack the obtained time series, and normalize it by the number of sites. The obtained time series reflect the average displacement over 1 day. We apply this procedure independently to a selection of 16 sites located in the area of the mergings (stations represented in blue in Supplementary Fig. 1) and to 10 sites away from the merging area in order to assess the level of noise (stations represented in red in Supplementary Fig. 1). We verify that sites located away from the SSE show spatially averaged daily velocity <0.2 mm day<sup>-1</sup> (Fig. 3b, red curve). On the contrary, sites located near the merging area show a clear increase in recorded slip displacement during the periods of enhanced tremor activity. In particular, the late-September merging phase is associated with a sharp (8x) increase in daily displacement (Fig. 3b, blue curve). This analysis demonstrates, independently from the inversion, that the 2013 sequence had several periods of coherent displacement acceleration, the largest one corresponding to the period of SSE merging around September 28 as found in our inversion.

## Data availability

The GPS time series and tremor catalog we used are respectively available on the UNAVCO (https://www.unavco.org) and PNSN websites (https://tremor.pnsn.org).

## Code availability

The pycs library we used for data processing and for the kinematic slip inversion is available upon request to J.-M.N. (email: nocquet@geoazur.unice.fr).

Received: 13 July 2019; Accepted: 27 February 2020; Published online: 01 May 2020

## References

1. Rogers, G. & Dragert, H. Episodic tremor and slip on the cascadia subduction zone: the chatter of silent slip. *Science* **300**, 1942–1943 (2003).

2. Obara, K. Nonvolcanic deep tremor associated with subduction in southwest Japan. *Science* **296**, 1679–1681 (2002).

3. Shelly, D. R., Beroza, G. C. & Ide, S. Non-volcanic tremor and low-frequency earthquake swarms. *Nature* **446**, 305 (2007).
4. Bartlow, N. M., Miyazaki, S., Bradley, A. M. & Segall, P. Space-time correlation of slip and tremor during the 2009 cascadia slow slip event. *Geophys. Res. Lett.* https://doi.org/10.1029/2011GL048714 (2011).

5. Wech, A. G. & Bartlow, N. M. Slip rate and tremor genesis in cascadia. *Geophys. Res. Lett.* **41**, 392–398 (2014).

6. Kao, H., Shan, S.-J., Dragert, H. & Rogers, G. Northern cascadia episodic tremor and slip: a decade of tremor observations from 1997 to 2007. *J. Geophys. Res. Solid Earth* https://doi.org/10.1029/2008JB006046 (2009).

7. Schmidt, D. & Gao, H. Source parameters and time-dependent slip distributions of slow slip events on the cascadia subduction zone from 1998 to 2008. *J. Geophys. Res. Solid Earth* https://doi.org/10.1029/2008JB006045 (2010).

8. Radiguet, M. et al. Spatial and temporal evolution of a long term slow slip event: the 2006 Guerrero slow slip event. *Geophys. J. Int.* **184**, 816–828 (2011).

9. Bartlow, N. M., Wallace, L. M., Beavan, R. J., Bannister, S. & Segall, P. Time-dependent modeling of slow slip events and associated seismicity and tremor at the Hikurangi subduction zone, New Zealand. *J. Geophys. Res. Solid Earth* **119**, 734–753 (2014).

10. Bletery, Q. et al. Characteristics of secondary slip fronts associated with slow earthquakes in cascadia. *Earth Planet. Sci. Lett.* **463**, 212–220 (2017).

11. Michel, S., Gualandi, A. & Avouac, J.-P. Interseismic coupling and slow slip events on the cascadia megathrust. *Pure Appl. Geophys.* **176**, 3867–3891 (2019).

12. Michel, S., Gualandi, A. & Avouac, J.-P. Similar scaling laws for earthquakes and cascadia slow-slip events. *Nature* **574**, 522–526 (2019).

13. Wech, A. G. & Creager, K. C. Automated detection and location of cascadia tremor. *Geophys. Res. Lett.* **35**, L20302 (2008).

14. Wech, A. G. Interactive tremor monitoring. *Seismol. Res. Lett.* **81**, 664–669 (2010).

15. Scholz, C. H. Earthquakes and friction laws. *Nature* **391**, 37 (1998).

16. Marone, C. Laboratory-derived friction laws and their application to seismic faulting. *Annu. Rev. Earth Planet. Sci.* **26**, 643–696 (1998).

17. Thomas, A. M., Nadeau, R. M. & Bürgmann, R. Tremor-tide correlations and near-lithostatic pore pressure on the deep San Andreas fault. *Nature* **462**, 1048 (2009).

18. Frank, W. B., Rousset, B., Lasserre, C. & Campillo, M. Revealing the cluster of slow transients behind a large slow slip event. *Sci. Adv.* **4**, eaat0661 (2018).

19. Thomas, A., Beeler, N., Bletery, Q., Burgmann, R. & Shelly, D. Using low-frequency earthquake families on the San Andreas fault as deep creepmeters. *J. Geophys. Res. Solid Earth* **123**, 457–475 (2018).

20. Dieterich, J. H. Modeling of rock friction: 2. Simulation of preseismic slip. *J. Geophys. Res. Solid Earth* **84**, 2169–2175 (1979).

21. Kaneko, Y. & Ampuero, J.-P. A mechanism for preseismic steady rupture fronts observed in laboratory experiments. *Geophys. Res. Lett.* **38**, L21307 (2011).

22. Fukuyama, E. et al. Spatiotemporal complexity of 2-d rupture nucleation process observed by direct monitoring during large-scale biaxial rock friction experiments. *Tectonophysics* **733**, 182–192 (2018).

23. Kato, A. et al. Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-oki earthquake. *Science* **335**, 705–708 (2012).

24. Radiguet, M. et al. Triggering of the 2014 Mw 7.3 Papanoa earthquake by a slow slip event in Guerrero, Mexico. *Nat. Geosci.* **9**, 829 (2016).

25. Voss, N. et al. Do slow slip events trigger large and great megathrust earthquakes. *Sci. Adv.* **4**, eaat8472 (2018).

26. Ruiz, S. et al. Intense foreshocks and a slow slip event preceded the 2014 Iquique Mw 8.1 earthquake. *Science* **345**, 1165–1169 (2014).

27. Kato, A. & Nakagawa, S. Multiple slow-slip events during a foreshock sequence of the 2014 Iquique, Chile Mw 8.1 earthquake. *Geophys. Res. Lett.* **41**, 5420–5427 (2014).

28. Herring, T. A. et al. Plate boundary observatory and related networks: GPS data analysis methods and geodetic products. *Rev. Geophys.* **54**, 759–808 (2016).

29. Hayes, G. P. et al. Slab2, a comprehensive subduction zone geometry model. *Science* **362**, 58–61 (2018).

30. DeMets, C., Gordon, R. G. & Argus, D. F. Geologically current plate motions. *Geophys. J. Int.* **181**, 1–80 (2010).

31. Meade, B. J. Algorithms for the calculation of exact displacements, strains, and stresses for triangular dislocation elements in a uniform elastic half space. *Comput. Geosci.* **33**, 1064–1075 (2007).

32. Stark, P. B. & Parker, R. L. Bounded-variable least-squares: an algorithm and applications. *Comput. Stat.* **10**, 129–129 (1995).

33. Nocquet, J.-M. Stochastic static fault slip inversion from geodetic data with non-negativity and bound constraints. *Geophys. J. Int.* **214**, 366–385 (2018).

34. Tarantola, A. *Inverse Problem Theory and Methods for Model Parameter Estimation*, Vol. 89 (Siam, 2005).

35. Nocquet, J.-M. et al. Motion of continental slivers and creeping subduction in the Northern Andes. *Nat. Geosci.* **7**, 287 (2014).

36. Vaca, S., Vallée, M., Nocquet, J.-M., Battaglia, J. & Régnier, M. Recurrent slow slip events as a barrier to the northward rupture propagation of the 2016 Pedernales earthquake (Central Ecuador). *Tectonophysics* **724**, 80–92 (2018).

37. Barbero, A. & Sra, S. Modular proximal optimization for multidimensional total-variation regularization. *J. Mach. Learn. Res.* **19**, 1–82 (2018).

## Acknowledgements

We built our fault geometry based on the slab2.0 model available online at the following address: https://www.sciencebase.gov/catalog/item/5aa1b00ee4b0b1c392e86467/. We thank A. Wech, K. Creager, and PNSN for letting their tremor catalog available (https://tremor.pnsn.org). We downloaded the GPS time series from the UNAVCO website (https://www.unavco.org). We used the Generic Mapping Tools GMT5 (https://gmt.soest.hawaii.edu/). We thank Jean-Paul Ampuero and Martin Vallée for helpful discussion. This work has been supported by the French government, through the UCAJEDI Investments in the Future project managed by the National Research Agency (ANR) ANR-15-IDEX-01 and by the ANR research grants E-POST, REMAKE, and S5.

## Author contributions

Q.B. had the original idea. J.-M.N. provided numerical tools (PYACS python library) and expertise. Both authors contributed to the analysis, interpretation and to the manuscript preparation.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41467-020-15494-4.

Correspondence and requests for materials should be addressed to Q.B.

Peer review information *Nature Communications* thanks Aitaro Kato and the other anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Reprints and permission information is available at http://www.nature.com/reprints

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_51.png" /></figure>

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020
## 4 Perspectives

Bien que l'étude des séismes lents est une piste intéressante en vue de la prédiction des séismes, nous sommes très loin d'une possible application opérationnelle. L'émergence de l'IA en sismologie apporte des possibilités nouvelles, peut-être plus prometteuses (voir Chapitre I.5). Le projet ERC EARLI, que je porte depuis Janvier 2021, comporte un volet très exploratoire sur les possibilités d'utilisation de l'IA à des fins prédictives. L'idée est d'exploiter une information de taille que nous avons à notre disposition que l'IA peut potentiellement utiliser efficacement : nous savons quand les séismes passés se sont produits. Cette information nous permet de définir la présence d'un signal précurseur dans les données selon l'occurrence postérieure d'un tremblement de terre. On peut théoriquement, à partir de là, entrainer un algorithme à reconnaître un précurseur dans le bruit sismique sans avoir besoin de l'identifier formellement (Figure III.2).

<figure><img src="image_52.png" /><figcaption>**Figure III.2** – Schéma illustratif du principe de la détection de précurseurs par apprentissage supervisé (voir Chapitre I.5.1). Les données d'entrainement sont composées de sismogrammes (3 composantes) de durées égales $t_w$. La moitié du jeu de données d'entrainement contient du bruit enregistré sur une station quelconque à un moment quelconque; le label correspondant à ces données est 0. L'autre moitié du jeu de données est composée de sismogrammes précédant un séisme d'un temps $t_p$; le label correspondant à ces données est 1. Une fois entrainé, l'algorithme prédira 0 ou 1 pour des sismogrammes sur lesquels il n'a pas été entrainé, 0 signifiant « rien à signaler », 1 signifiant « séisme attendu dans $t_p$ ».</figcaption></figure>

Quelque soit l'origine des signaux précurseurs, on s'attend à ce qu'ils dépendent, au premier
ordre, de la distance $d_{i,k}$ entre la source $S^i$ et la station $ST^k$ et de la magnitude $M_w^i$. Comme on connaît la localisation des séismes passés, on peut facilement calculer $d_{i,k}$, et l'utiliser comme label avec $M_w^i$ (Figure III.3).

<figure><img src="image_53.png" /><figcaption>**Figure III.3** – Schéma illustratif de la modification du protocole d'entrainement par rapport à la version simplifiée (Figure III.2). Le label binaire est augmenté avec les informations sur la magnitude $M_w^i$ de l'évènement et sa distance $d_{i,k}$ à la station $ST^k$. Lors de la phase prédictive, l'algorithme pourrait alors procurer des informations sur la magnitude et la distance de l'évènement à venir, permettant ainsi de prédire la localisation (en croisant les informations données par différentes stations) et la magnitude.</figcaption></figure>

L'étape suivante serait alors d'inclure $t_p$ dans le label, afin d'avoir également accès au temps avant l'évènement lors de la phase prédictive. L'approche proposée ici est clairement explo- ratoire et repose sur (1) l'hypothèse de l'existence des précurseurs et (2) sur le fait que l'on dispose de suffisamment de données pour extraire ces signaux. Notez, toutefois, qu'une grande force de cette approche est de pouvoir, en théorie, fonctionner quelque soit l'origine des signaux précurseurs. Qu'ils soient générés par une phase de glissement lent précurserice, par une coales- cence de séismes lents, par un changement des paramètres physiques du milieu ou quelqu'autre phénomène physique qu'il soit, nous n'avons besoin ici de ne faire aucune hypothèse sur l'origine des précurseurs.
# Chapitre IV

## Comment anticiper les grands séismes ?

### 1 L’alerte sismique/tsunami, une course contre le temps

Nous avons vu dans les chapitres précédents que, si quelques pistes existent, la prévision opérationnelle des séismes est difficilement envisageable à moyen terme. Néanmoins, il n’est pas forcément nécessaire de prédire le futur pour anticiper les catastrophes. Les systèmes d’alerte sismique détectent et caractérisent rapidement les séismes en cours en temps réel dans le but d’alerter les populations alentours de l’imminence de secousses. Ils utilisent l’information contenue dans les premières secondes d’enregistrement de l’évènement, qui sont généralement d’assez faible amplitude, pour estimer de potentielles secousses de plus forte amplitude. Les messages d’alerte se transmettant plus vite que les ondes sismiques ne se propagent, ils peuvent atteindre les populations concernées avant les premières secousses, offrant un temps précieux aux personnes et aux systèmes automatisés pour prendre des mesures afin de réduire les pertes humaines et matérielles. Les mesures vont de simples procédures invitant la population à rejoindre un lieu sûr à des procédures automatisées complexes comme interrompre le décollage et l’atterrissage des avions [Wu et al., 2013].

Des systèmes d’alerte opèrent dans de nombreuses régions du monde exposées aux tremblements de terre. Le Mexique a un système d’alerte opérationnel depuis 1991 [Espinosa-Aranda et al., 2009]. L’Agence Météorolgique du Japon (JMA) opère un système depuis 2007 [Hoshiba et al., 2011]. Comme le système japonais, le système ShakeAlert développé pour la Côte Ouest des États-Unis [Kohler et al., 2018] estime la localisation et la magnitude des séismes. Ces estimations sont alors utilisées dans des équations de prédiction des mouvements du sol pour estimer l’amplitude des secousses ainsi que le risque de tsunami; une alerte est lancée – ou non – sur la base de ces dernières estimations. Le système JMA a lancé des centaines d’alertes, dont certaines ont concerné des millions de personnes, comme pour le séisme de Tohoku-Oki (Mw 9.0) en 2011 [Fujinawa and Noda, 2013].
Les systèmes d'alerte ont été conçus dans l'hypothèse que la taille et l'étendue de la rupture pouvaient être connues immédiatement [Heaton, 1985; Meier, 2017] ou que la magnitude finale d'un séisme pouvait être déterminée à partir des premières secondes d'enregistrement des ondes P [Allen, 2006; Zollo et al., 2006]. Or, des études indiquent que la magnitude finale d'un séisme n'est pas déterminée au début de la rupture [Kanamori, 2005; Rydelek and Horiuchi, 2006; Rydelek et al., 2007; Yamamoto et al., 2008]. Des arguments théoriques vont dans ce sens [Sato and Kanamori, 1999] ainsi que des études observationnelles qui suggèrent que les premières secondes des sismogrammes des grands et des petits séismes sont indissociables [Mori and Kanamori, 1996; Ellsworth and Beroza, 1998; Kilb and Gomberg, 1999; Uchide and Ide, 2010; Meier et al., 2017].

<figure><img src="image_54.png" /><figcaption>**Figure IV.1** – Estimation de la magnitude par le système JMA (bleu) en fonction du temps après l'origine des séismes de (A) Tohoku-OKi (M<sub>w</sub> 9.0, 2011), (B) Kumamoto (M<sub>w</sub> 7.0, 2016) et (C) Iwate-Miyagi Nairiku (M<sub>w</sub> 6.9, 2008) [Minson et al., 2018]. Les courbes noires donnent les fonctions sources des séismes ( évolution réelle de la magnitude en fonction du temps).</figcaption></figure>

L'observation de la performance des systèmes d'alerte opérationnels va également dans ce sens (Figure IV.1). L'estimation de la magnitude de séismes récents, tels que les 3 présentés dans la Figure IV.1, suggère que l'estimation suit l'évolution du moment relâché par le séisme avec un délai dû à l'attente de l'arrivée des premières ondes sismiques. On note également que le système JMA sature au delà de M<sub>w</sub> = 8.1. L'étude présentée ci-après a pour objectif de réduire le délai de latence (dû à l'attente de l'arrivée des ondes sismiques) et de résoudre le problème
de saturation, problématique en particulier pour l'estimation du tsunami, comme dans le cas
du séisme de Tohoku-Oki.

2 Estimation instantanée de la croissance des très grands séismes

L'étude qui suit est venue de l'idée de tester l'approche décrite dans les perspectives du Chapitre précédent (section III.4) sur des signaux intermédiaires, non précurseurs puisqu'ils sont générés par le séisme, mais précédant les ondes sismiques puisqu'ils présentent la particularité de se propager à la vitesse de la lumière : les PEGS (voir Chapitre I.4.3). Nous avons ainsi conçu un algorithme d'IA, PEGSnet, capable d'estimer la magnitude et la localisation de très grands séismes à partir des PEGS enregistrés sur un réseau de stations sismiques. Cet algorithme, entrainé sur un très grand nombre de données synthétiques augmentées de bruit empirique, s'avère capable d'estimer la magnitude des très grands séismes de manière instantanée – i.e. sans le délai dû à l'attente des ondes sismiques – à partir du moment où cette magnitude atteint Mw = 8.3. S'il n'est pas performant pour des évènements de magnitude inférieure à 8.3, PEGSnet résout le problème de saturation des très grandes magnitudes. Ses performances sont d'autant plus impressionnantes que PEGSnet n'utilise que des données précédant les premières ondes sismiques, i.e. exclusivement des formes d'onde communément traitées comme du bruit.

Cette étude a été réalisée à Géoazur dans le cadre du postdoctorat d'Andrea Licciardi,
en collaboration avec Bertrand Rouet-Leduc, Jean-Paul Ampuero et Kévin Juhel. Elle a été
publiée en 2022 dans le journal Nature [Licciardi et al., 2022].
# Instantaneous tracking of earthquake growth with elastogravity signals

https://doi.org/10.1038/s41586-022-04672-7

Received: 14 June 2021

Accepted: 22 March 2022

Published online: 11 May 2022

Open access

<figure><img src="image_55.png" /></figure>

Check for updates

Andrea Licciardi<fnref n="1" />,²,<fnref n="3" />, Quentin Bletery<fnref n="1" />, Bertrand Rouet-Leduc²,<fnref n="3" />, Jean-Paul Ampuero<fnref n="1" /> & Kévin Juhel<fnref n="1" />,⁴

Rapid and reliable estimation of large earthquake magnitude (above 8) is key to mitigating the risks associated with strong shaking and tsunamis<fnref n="1" />. Standard early warning systems based on seismic waves fail to rapidly estimate the size of such large earthquakes²⁻<fnref n="5" />. Geodesy-based approaches provide better estimations, but are also subject to large uncertainties and latency associated with the slowness of seismic waves. Recently discovered speed-of-light prompt elastogravity signals (PEGS) have raised hopes that these limitations may be overcome⁶,⁷, but have not been tested for operational early warning. Here we show that PEGS can be used in real time to track earthquake growth instantaneously after the event reaches a certain magnitude. We develop a deep learning model that leverages the information carried by PEGS recorded by regional broadband seismometers in Japan before the arrival of seismic waves. After training on a database of synthetic waveforms augmented with empirical noise, we show that the algorithm can instantaneously track an earthquake source time function on real data. Our model unlocks ‘true real-time’ access to the rupture evolution of large earthquakes using a portion of seismograms that is routinely treated as noise, and can be immediately transformative for tsunami early warning.

The sudden displacement of rock mass induced by an earthquake generates density variations that, in turn, modify the Earth's gravity field. The signal associated with these transient gravity perturbations propagates at the speed of light, much faster than the fastest elastic waves (P-waves)⁸⁻¹⁰. Recent theoretical studies have shown the potential for earthquake early warning systems (EEWS) that are based on the gravity signals that would be measured by a future generation of gravity-gradient sensors¹¹,¹², yet to be developed. In practice, existing inertial sensors (for example, seismometers) measure a combination of the direct gravity perturbations and their induced elastic response, named prompt elastogravity signals (PEGS)¹³,¹⁴. PEGS detection on real data is difficult for two reasons. First, the amplitude of the direct gravity perturbations is very small. Second, the induced elastic response tends to cancel out the gravity effects on seismometer recordings, especially in the early portion of the signal. The combination of these effects results in detectability limited to a time window preceding the P-wave arrival, which depends on epicentral distance (between a few seconds to a few tens of seconds), where PEGS reach their maximum amplitudes (a few nm s⁻² at most)¹⁴,¹⁵. Nevertheless, PEGS could prove beneficial for EEWS. First, they travel at the speed of light and might provide extra time for alert. Second, PEGS do not saturate, as opposed to P-waves recorded by near-field seismometers that may clip during large earthquakes. Finally, given the wavelength of the signal and the smoothness of the generated wavefield¹⁴, the spatial complexity of the rupture does not substantially affect PEGS amplitudes¹⁵. For this reason, PEGS dependence on earthquake magnitude, focal mechanism

and source time function (STF) has the potential to improve early characterization of earthquake size and source parameters, under a simple point-source approximation¹⁵,¹⁶. In this work we show that PEGS can efficiently be used to improve operational EEWS.

Given the expected level of background noise, previous works have suggested a limit for PEGS detectability to earthquakes with magnitude (*M*<sub>w</sub> above 8 (refs. ¹⁵,¹⁶). The occurrence of *M*<sub>w</sub> > 8 earthquakes poses a difficult challenge for conventional EEWS. On the one hand, subduction megathrust earthquakes require accurate and fast estimates of final magnitude to mitigate the risk associated with strong shaking and to forecast the potential size of tsunami waves¹⁷⁻¹⁹. On the other hand, EEWS based on point-source algorithms that rely on the first few seconds of P-waves tend to produce saturated magnitude estimation for such large earthquakes. One reason is that instruments saturate for very large events. Another more fundamental reason is that the early portion of seismograms simply does not contain enough information to distinguish between small and large earthquakes (which have longer duration) at the very early stage of rupture<fnref n="5" />. An example of this paradigm is the performance of the EEWS of the Japan Meteorological Agency (JMA) during the 2011 *M*<sub>w</sub> = 9.0 Tohoku-Oki earthquake, which underestimated the final *M*<sub>w</sub> of the event to around 8.1 (ref. ²⁰). Although the deterministic nature of earthquake rupture is still debated, a growing amount of evidence suggests that earthquake ruptures are not deterministic²¹. Therefore, EEWS should be designed to track the moment release as the rupture unfolds instead of forecast the final earthquake magnitude<fnref n="3" />,²². Over the last decade, finite-fault EEWS based on global navigation satellite

<footnote n="1">¹Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Géoazur, Sophia Antipolis, France. ²Disaster Prevention Research Institute, Kyoto University, Kyoto, Japan.</footnote>

<footnote n="3">³Los Alamos National Laboratory, Geophysics Group, Los Alamos, NM, USA. ⁴Laboratoire de Planétologie et Géodynamique, UMR 6112, Nantes University, CNRS, Nantes, France.</footnote>

<footnote n="5">⁵e-mail: andrea.licciardi@geoazur.unice.fr</footnote>
<figure><img src="image_56.png" /><figcaption>**Fig. 1** | Experimental set-up and input data examples from the synthetic database. a, Locations of the stations in our network are indicated by light red symbols. Triangles are stations belonging to the F-Net seismic network and squares are stations from various networks obtained through the IRIS data centre. Dark green lines are made of 1,400 discrete locations of the synthetic earthquakes considered in our training database (see text). The orange star indicates the source location for one random synthetic example ($M_w$ = 8.9) extracted from the database whose moment rate (blue) and $M_w$ evolution (grey)</figcaption></figure>

<figure><img src="image_57.png" /></figure>

are shown in the inset. b, The corresponding vertical waveforms at stations MDJ and YSI, where the red trace is the noise-free synthetic PEGS and the black trace is the same, with empirical noise added. Only the pre-P-wave time window is shown. c, Synthetic vertical waveforms plus noise associated with the selected event for each station in the seismic network. The grey dashed line indicates the P-wave arrival time. Stations are sorted by their longitude. The positions of MDJ and YSI are indicated by black horizontal dashed lines.

system (GNSS) data have emerged as a new tool with which to overcome the magnitude saturation problem<sup>23–28</sup>. Nevertheless, subjective choices required in GNSS data selection and/or preprocessing<sup>24,25</sup> may result in large uncertainties. Moreover, the fast responses achieved for megathrust earthquakes<sup>26</sup> have recently been questioned and attributed not to the predictive power of GNSS data—which would enable the estimation of an earthquake’s final magnitude before the rupture is over—but to prior constraints and regularization-induced artefacts<sup>22</sup>. A deep learning model based on GNNS data has recently been proposed to overcome these limitations<sup>29</sup>. Although it proved promising, as for other finite-fault approaches, it requires a priori assumptions on slip distribution. Finally, all existing EEWS suffer from unavoidable latency, owing to the speed at which the information is carried by P-waves, and therefore produce a time-shifted version of the earthquake STF.

In this context, we show that a convolutional neural network (CNN)<sup>30</sup> approach can leverage the information carried by PEGS at the speed of light to overcome these limitations for large earthquakes. Successful applications of deep learning in seismology have provided new tools for pushing the detection limit of small seismic signals<sup>31,32</sup> and for the characterization of earthquake source parameters (magnitude and location)<sup>33–35</sup> with EEWS applications<sup>29,36,37</sup>. Here we present a deep learning model, PEGSNet, trained to estimate earthquake location and track the time-dependent magnitude, $M_w(t)$, from PEGS data before P-wave arrivals.

## The training database

We focus on the Japanese subduction zone because: (1) it has experienced one of the largest and most destructive subduction events (the 2011 $M_w$ = 9.0 Tohoku-Oki earthquake), and (2) we can rely on a large and dense network of good-quality broadband seismic stations (Fig. 1a).

Because of the paucity of PEGS observations<sup>16</sup>, we train PEGSNet on a database of synthetic PEGS waveforms generated with a normal-mode summation code<sup>14</sup> and add real noise recorded at each station in the seismic network. The training database is made of 500,000 synthetic earthquake sources distributed along the Japanese megathrust, with location, strike and dip angles following the Slab2.0 model<sup>38</sup>. We draw random magnitude and rake angles following uniform (from 5.5 to 10) and Gaussian ($\mu$ = 90°, $\sigma$ = 10°) distributions, respectively. In this work, we consider only the most simplistic rupture spatial descriptions: point sources. The difference in PEGS predictions associated with the source finiteness was shown to be within data uncertainties during the 2011 Tohoku earthquake<sup>15</sup>. The validity of this assumption is further corroborated by the detection of PEGS generated by multiple events considering a simple point-source model<sup>16</sup>. However, simple triangular STFs might not suffice for modelling PEGS for realistic earthquake ruptures<sup>16</sup>. For this reason, we generate random STFs using a model designed to mimic empirical laws and statistical observations<sup>5</sup> (Fig. 1a, inset, and Extended Data Fig. 1a) to produce
<figure><img src="image_58.png" /><figcaption>**Fig. 2 | Tracking M<sub>w</sub>(t) from PEGS data.** Results of predictions on the test set. **a,** Predictions accuracy as a function of time and M<sub>w</sub>. For each pixel in the image, accuracy is calculated as the number of successful predictions divided by the total number of samples. A prediction is defined as successful if the distance with its corresponding label is within ±0.4 magnitude units. Dashed lines indicate the average values of the true M<sub>w</sub>(t) for events with different final</figcaption></figure>

magnitudes. **b,** As in **a,** but showing the average residuals ($|M_{w[true]} - M_{w[pred]}|$) for each pixel. **c,** Probability density of M<sub>w</sub> predictions for all the events in the test set with true final M<sub>w</sub> of 9 ± 0.05. The solid red line is the mode of the distribution. The red dashed lines bound the 25th–95th percentile range. The blue lines are the median (solid) and the 5th–95th percentile range (dashed) of the labels.

three-component (east, north, vertical) waveforms representing the source characteristics of all expected large megathrust earthquakes. In this STF model the smaller and larger earthquakes start in the same way. We then add empirical noise recorded at each station (Fig. 1b) to the generated waveforms. Both synthetics and noise seismograms are decimated to 1 Hz and bandpass-filtered between 2 mHz (high-pass Butterworth, two poles, causal) and 30 mHz (low-pass Butterworth, six poles, causal), to enhance long-period PEGS amplitudes by suppressing higher-frequency noise<sup>7</sup>. Traces are clipped to a threshold value of ±10 nm s<sup>-2</sup> (the maximum expected amplitude for PEGS in our database) and scaled by the same value to facilitate convergence of the optimizer, and at the same time to preserve information about PEGS radiation patterns across the seismic network. Finally, we store each example as 700-s-long traces centred at the origin time of the earthquake (Fig. 1c), with its three labels corresponding to the latitude (φ), longitude (λ) and the time-integrated STF of the event converted to the time-dependent magnitude M<sub>w</sub>(t).

## Building PEGSNet, labelling and training

PEGSNet is a deep CNN that combines convolutional layers and fully connected layers in sequence (Extended Data Fig. 2; details in Methods). There exists a wealth of literature on deep learning architectures, but our goal here is to show the feasibility of PEGS-based EEWS using deep learning, and we therefore use a simple and classic Le-NET-like CNN architecture<sup>39</sup>. Similar architectures have been used for magnitude and location estimation using P and S-waves with a single-station approach to increase the number of training examples<sup>33,34</sup>. Here we take advantage of our large synthetic database and adopt a multi-station approach.

CNNs work best with image-like inputs, and we arrange our data as an image, sorting stations by their longitude (Fig. 1c).

Our aim is to track the time-dependent magnitude, M<sub>w</sub>(t), which corresponds to the time-integrated STF at a given time. For PEGSNet to learn M<sub>w</sub>(t), we design a specific strategy: we randomly select the starting time (T<sub>1</sub>) of the analysed data during training. For each input example, we extract 315 s of data from T<sub>1</sub> to T<sub>2</sub> = T<sub>1</sub> + 315 s and assign the value of M<sub>w</sub>(T<sub>2</sub>) as the corresponding label. This allows the model to learn patterns in the data as the STF evolves with time.

We exclusively use the information carried by PEGS by setting to zero the amplitude of the seismograms for t ≥ T<sub>p</sub>, where T<sub>p</sub> is the P-wave arrival time at each station for a given event. Therefore, P-wave arrival times are assumed to be known—in practice, any existing triggering algorithms (whether based on deep learning or not) could be used to obtain them before feeding the data to the model. The final input data for each example consists of a three-channel image (one for each seismic component), and the output layer of PEGSNet predicts three real values for M<sub>w</sub>(T<sub>2</sub>), φ and λ.

The synthetic database described above is randomly split into 70% training, 20% validation and 10% test sets, with a given synthetic STF and a given empirical noise measurement going to only one of the datasets. Training is performed by minimizing the Huber loss between the labels (the true values of M<sub>w</sub>(t), φ and λ are known) and the predicted values. The model with the lowest error on the validation set is chosen as the final model.

## Instantaneous tracking of M<sub>w</sub>

Figure 2 summarizes the results on the synthetic test set for M<sub>w</sub>(t), and those for φ and λ are displayed in Extended Data Fig. 3. Results
<figure><img src="image_59.png" /><figcaption>**Fig. 3** | Moment tracking of the 2011 M<sub>w</sub> 9.0 Tohoku-Oki earthquake. **a,** Comparison of PEGSNet M<sub>w</sub>(t) predictions (blue) with the source time function (STF) of the Tohoku-Oki earthquake<sup>40</sup> (orange) and with the results of three existing EEWS (JMA<sup>20</sup>, FinDer2<sup>41</sup> and BEFORES<sup>24</sup>). Dashed orange lines indicate</figcaption></figure>

±0.3 magnitude units. **b,** Difference between the predicted M<sub>w</sub> of various algorithms (PEGSNet in blue) and the 'true' STF as shown in **a**. The solid orange line indicates a difference of zero, the dashed orange line a difference of -0.3 M<sub>w</sub> units (underestimation).

for location rely on the P-wave arrival times (which are assumed to be known) and show errors around 25 to 30 km starting at about 50 s after origin. In the following, we will exclusively discuss the results for M<sub>w</sub>(t). These results represent the main contribution of this work as they are obtained using only the pre P-wave portion of the seismograms, which is routinely treated as pure noise. For each example in the test set, we simulate a real-time scenario by parsing the data with a running time window of 315 s. Each data window is fed to PEGSNet, which makes an M<sub>w</sub> estimate for the end of the window, and the window is shifted in steps of 1 s as PEGSNet attempts to progressively reconstruct the STF (Extended Data Fig. 4). We define a successful prediction if the estimated M<sub>w</sub>(t) lies within ±0.4 magnitude units from the true value.

PEGSNet can track the moment released by earthquakes with final M<sub>w</sub> above 8.6 with good accuracy (above 90%, Fig. 2a) and low errors (M<sub>w</sub> error below 0.25, Fig. 2b), starting at about 40 s after origin time. For earthquakes with final M<sub>w</sub> between 8.2 and 8.6, early tracking of the moment release is more difficult and only the final M<sub>w</sub> can be estimated (with accuracy between 60% and 70% and errors above 0.25) after 150 s from origin time. Any predicted M<sub>w</sub> below 8.2 is poorly constrained by the data and we place a conservative lower limit on PEGSNet sensitivity to M<sub>w</sub> at 8.3. However, this limit depends on the noise amplitudes across the seismic network. Under favourable noise conditions (standard deviation of the noise for the whole seismic network below around 0.5 nm s<sup>-2</sup>), it can be reduced to about 7.9 or 8.0 (Extended Data Fig. 5).

To examine the time-dependent performance of PEGSNet more closely, Fig. 2c shows all the predictions (on the test set) associated with events with a final M<sub>w</sub> of 9.0 ± 0.05. In the first 30–40 s after origin, predictions are strongly underestimated. This arises from a combination of two effects: the lack of sensitivity below 8.3 (in the first 30–40 s after origin the current magnitude is below 8.3 in our STF database) and the fact that, at stations located close to the source, P-waves can mask PEGS, owing to their small time separation. A test performed on noise-free synthetics shows that at least 15 s after origin are needed by the model to reliably predict M<sub>w</sub> (Extended Data Fig. 6). Starting at about 40 s (when the true M<sub>w</sub> is generally above 8.2 in our STF database), PEGSNet can track the evolution of the moment release instantaneously, that is, without time shift between estimated and true M<sub>w</sub>(t), as indicated by the mode of the distribution of the predictions. Our results show that PEGSNet can exploit the key features of PEGS: (1) PEGS sensitivity to magnitude for large earthquakes allows the model to distinguish, for example, between an M<sub>w</sub> = 8 and an M<sub>w</sub> = 9 earthquake and (2) the information about M<sub>w</sub> is effectively propagated at the speed of light,

which results in instantaneous tracking of the moment release, from about 40 s into the rupture until its end. Although P-wave triggering is needed at each station, the instantaneous information about the source carried by PEGS is readily available in the data that precede P-wave arrivals. This allows PEGSNet to estimate M<sub>w</sub>(t) with zero delay once the magnitude exceeds approximately 8.3.

## Playback of the Tohoku-Oki earthquake

We tested PEGSNet on real data from the 2011 M<sub>w</sub> 9.0 Tohoku-Oki earthquake. The raw seismograms are processed as described in Methods and fed to PEGSNet. We simulated a data playback scenario and made an M<sub>w</sub> prediction at each second starting at the earthquake origin time and using only the previous 315 s of data (details in Methods). Figure 3 shows the results of this retrospective analysis compared with a 'true' STF calibrated against various different data types<sup>40</sup>, the historical JMA EEWS performance<sup>20</sup>, and the best-performing finite-fault EEWS based on seismic (FinDer2)<sup>41</sup> and GNSS data (BEFORES)<sup>24</sup>. Further comparisons with additional algorithms are shown in Extended Data Fig. 7. After around 55 s, PEGSNet is always closer to the 'true' STF than the other algorithms (Fig. 3a). The predicted M<sub>w</sub> between 55 and 100 s is underestimated by about 0.3 magnitude units compared to the ground truth. PEGSNet reaches a correct prediction of the final M<sub>w</sub> after about 120 s, when the rupture is almost over. Despite the slight underprediction, the predicted M<sub>w</sub>(t) values between 55 and 100 s show a clear increasing trend and suggest that the rupture is still in progress. Therefore, we argue that PEGS can be used in real time to track the evolving rupture as it unfolds. Consistently with our synthetic tests (Fig. 2), our predictions after about 55 s (corresponding to the time the earthquake reaches M<sub>w</sub> ≈ 8.3) provide a snapshot of the evolution of the earthquake rupture at that exact time because PEGS carry information about the source at the speed of light. This is further corroborated by a test on synthetic data generated using the 'true' Tohoku-Oki earthquake STF and noise seismograms extracted from the test set (Extended Data Fig. 8).

Notably, PEGSNet does not suffer from the magnitude saturation issue that instead affects the seismic-based EEWS. Both in terms of latency and accuracy, the performance of the GNSS-based algorithm BEFORES is the closest to PEGSNet. However, the response of PEGSNet is faster: for an arbitrary threshold at M<sub>w</sub> = 8.0, PEGSNet would issue an alert 53 s after origin, BEFORES about 8 s later. PEGSNet is more accurate: the difference with the 'true' STF is always equal to or less than 0.3 M<sub>w</sub> units starting at 55 s, whereas BEFORES reaches that level
of precision only after 100 s from origin (Fig. 3b). More importantly, in
contrast to GNSS-based EEWS, PEGSNet does not require prior assump-
tions on the data, because it is trained to learn features from the entire
network of stations regardless of their data quality. In addition, no-slip
distributions are assumed or modelled in our approach as PEGS are only
sensitive to smooth variations of the released moment and, therefore,
less likely to be influenced by the spatial complexity of the rupture.

To further assess the robustness of PEGSNet predictions on the
Tohoku data, we performed a test in which we keep the P-wave
arrival information but substitute the actual recorded waveforms
in the pre-P-wave time window with noise. The resulting prediction
of M<sub>w</sub>(t) never exceeds the lower sensitivity limit of PEGSNet (that
is, 8.3) and remains constant at about 6.5, which provides a baseline
value for noise (Extended Data Fig. 9). We also tested PEGSNet on all
the subduction earthquakes (dip-slip mechanism within 40 km from
the megathrust) with M<sub>w</sub> ≥ 7 that have occurred since January 2003,
without considering aftershocks (Methods and Extended Data Fig. 10).
For M<sub>w</sub> < 8 earthquakes, PEGSNet predictions converge toward the
noise baseline, confirming that these events are essentially not dis-
tinguishable from noise.

Implications for early warning

We have demonstrated instantaneous tracking of moment release for
large earthquakes (Fig. 3). Our results promote PEGS as a new class
of observables, easily accessible from the recordings of currently
deployed broadband seismometers worldwide, for practical applica-
tion in early warning systems that are currently limited by the speed of
P-waves. In the context of EEWS, PEGSNet can complement any
existing algorithm (either seismic- or GNSS-based), to improve M<sub>w</sub>
latency estimation and accuracy for M<sub>w</sub> > 8.3 megathrust earthquakes.
For example, PEGSNet could be combined with a recent deep learn-
ing model based on GNSS data²⁹ to eliminate intrinsic latency due to
P-wave speed.

At the same time, PEGSNet can immediately prove critical for tsunami early warning for which *M*<sub>w</sub> estimation within a few minutes is vital. Continuous updates of current *M*<sub>w</sub> can be fed to predictive models of tsunami waves height, helping mitigate the associated risk.

Our results suggest that PEGS can play a key part in the early charac-
terization of rupture evolution for large earthquakes. For such events,
PEGS data represent a new and independent source of information to
constrain the magnitude in real time.

PEGSNet requires only a few modifications to be implemented in
real time. Once trained, PEGSNet predictions are quasi-instantaneous,
although some latency could be introduced by the preprocessing step.
Although tailored here to earthquakes from the Japanese subduction
zone, PEGSNet can be easily adapted to other seismic networks and
regions and source mechanisms. In particular, PEGSNet portability to
other regions only requires the availability of noise recordings for the
seismic network of interest.

Online content

Any methods, additional references, Nature Research reporting sum-
aries, source data, extended data, supplementary information,
acknowledgements, peer review information; details of author contri-
butions and competing interests; and statements of data and code avail-
ability are available at https://doi.org/10.1038/s41586-022-04672-7.

1. Allen, R. M. & Melgar, D. Earthquake early warning: advances, scientific challenges, and societal needs. *Annu. Rev. Earth. Planet. Sci.* **47**, 361–388 (2019).

2. Wald, D. J. Practical limitations of earthquake early warning. *Earthq. Spectra* **36**, 1412–1447 (2020).

3. Minson, S. E., Meier, M. A., Baltay, A. S., Hanks, T. C. & Cochran, E. S. The limits of earthquake early warning: timeliness of ground motion estimates. *Sci. Adv.* **4**, eaaq0504 (2018).

4. Meier, M. A. How “good” are real-time ground motion predictions from earthquake early warning systems? *J. Geophys. Res. Solid Earth* **122**, 5561–5577 (2017).

5. Meier, M. A., Ampuero, J.-P. & Heaton, T. H. The hidden simplicity of subduction megathrust earthquakes. *Science* **357**, 1277–1281 (2017).

6. Montagner, J.-P. et al. Prompt gravity signal induced by the 2011 Tohoku-Oki earthquake. *Nat. Commun.* **7**, 13349 (2016).

7. Vallée, M. et al. Observations and modeling of the elastogravity signals preceding direct seismic waves. *Science* **358**, 1164–1168 (2017).

8. Dahlen, F. & Tromp, J. *Theoretical Global Seismology* (Princeton Univ. Press, 1998).

9. Harms, J. et al. Transient gravity perturbations induced by earthquake rupture. *Geophys. J. Int.* **201**, 1416–1425 (2015).

10. Harms, J. Transient gravity perturbations from a double-couple in a homogeneous half-space. *Geophys. J. Int.* **205**, 1153–1164 (2016).

11. Juhel, K. et al. Earthquake early warning using future generation gravity strainmeters. *J. Geophys. Res. Solid Earth* **123**, 10889–10902 (2018).

12. Shimoda, T., Juhel, K., Ampuero, J.-P., Montagner, J.-P. & Barsuglia, M. Early earthquake detection capabilities of different types of future-generation gravity gradiometers. *Geophys. J. Int.* **224**, 533–542 (2021).

13. Heaton, T. H. Correspondence: response of a gravimeter to an instantaneous step in gravity. *Nat. Commun.* **8**, 966 (2017).

14. Juhel, K. et al. Normal mode simulation of prompt elasto-gravity signals induced by an earthquake rupture. *Geophys. J. Int.* **216**, 935–947 (2019).

15. Zhang, S., Wang, R., Dahm, T., Zhou, S. & Heimann, S. Prompt elasto-gravity signals (PEGS) and their potential use in modern seismology. *Earth Planet. Sci. Lett.* **536**, 116150 (2020).

16. Vallée, M. & Juhel, K. Multiple observations of the prompt elastogravity signals heralding direct seismic waves. *J. Geophys. Res. Solid Earth* **124**, 2970–2989 (2019).

17. Melgar, D. & Bock, Y. Near-field tsunami models with rapid earthquake source inversions from land- and ocean-based observations: the potential for forecast and warning. *J. Geophys. Res. Solid Earth* **118**, 5939–5955 (2013).

18. Melgar, D. et al. Local tsunami warnings: perspectives from recent large events. *Geophys. Res. Lett.* **43**, 1109–1117 (2016).

19. Setiyono, U., Gusman, A. R., Satake, K. & Fujii, Y. Pre-computed tsunami inundation database and forecast simulation in Pelabuhan Ratu, Indonesia. *Pure Appl. Geophys.* **174**, 3219–3235 (2017).

20. Hoshiba, M. & Ozaki, T. Earthquake early warning and tsunami warning of the Japan Meteorological Agency, and their performance in the 2011 off the Pacific Coast of Tohoku Earthquake (M<sub>w</sub> 9.0). In *Early Warning for Geological Disasters* (eds Wenzel, F. & Zschau, J.) 1–28 (Springer, 2014).

21. Meier, M.-A., Ampuero, J.-P., Cochran, E. & Page, M. Apparent earthquake rupture predictability. *Geophys. J. Int.* **225**, 657–663 (2020).

22. McGuire, J. J., Minson, S. E., Murray, J. R. & Brooks, B. A. The role of geodetic algorithms for earthquake early warning in Cascadia. *Geophys. Res. Lett.* **48**, e2020GL092324 (2021).

23. Grapenthin, R., Johanson, I. A. & Allen, R. M. Operational real-time GPS-enhanced earthquake early warning. *J. Geophys. Res. Solid Earth* **119**, 7944–7965 (2014).

24. Minson, S. E., Murray, J. R., Langbein, J. O. & Gomberg, J. S. Real-time inversions for finite fault slip models and rupture geometry based on high-rate GPS data. *J. Geophys. Res. Solid Earth* **119**, 3201–3231 (2014).

25. Melgar, D. et al. Earthquake magnitude calculation without saturation from the scaling of peak ground displacement. *Geophys. Res. Lett.* **42**, 5197–5205 (2015).

26. Ruhl, C. J., Melgar, D., Grapenthin, R. & Allen, R. M. The value of real-time GNSS to earthquake early warning. *Geophys. Res. Lett.* **44**, 8311–8319 (2017).

27. Crowell, B. W. et al. G-FAST earthquake early warning potential for great earthquakes in Chile. *Seismol. Res. Lett.* **89**, 542–556 (2018).

28. Ruhl, C. J., Melgar, D., Chung, A. I., Grapenthin, R. & Allen, R. M. Quantifying the value of real-time geodetic constraints for earthquake early warning using a global seismic and geodetic data set. *J. Geophys. Res. Solid Earth* **124**, 3819–3837 (2019).

29. Lin, J. T., Melgar, D., Thomas, A. M. & Searcy, J. Early warning for great earthquakes from characterization of crustal deformation patterns with deep learning. *J. Geophys. Res. Solid Earth* **126**, e2021JB022703 (2021).

30. LeCun, Y., Bengio, Y. & Hinton, G. Deep learning. *Nature* **521**, 436–444 (2015).

31. Ross, Z. E., Meier, M.-A. & Hauksson, E. P wave arrival picking and first-motion polarity determination with deep learning. *J. Geophys. Res. Solid Earth* **123**, 5120–5129 (2018).

32. Rouet-Leduc, B., Hulbert, C., McBrearty, I. W. & Johnson, P. A. Probing slow earthquakes with deep learning. *Geophys. Res. Lett.* **47**, e2019GL085870 (2020).

33. Perol, T., Gharbi, M. & Denolle, M. Convolutional neural network for earthquake detection and location. *Sci. Adv.* **4**, e1700578 (2018).

34. Mousavi, S. M. & Beroza, G. C. A machine-learning approach for earthquake magnitude estimation. *Geophys. Res. Lett.* **47**, e2019GL085976 (2020).

35. van den Ende, M. P. A. & Ampuero, J.-P. Automated seismic source characterisation using deep graph neural networks. *Geophys. Res. Lett.* **47**, e2020GL088690 (2020).

36. Lomax, A., Michelini, A. & Jozinović, D. An investigation of rapid earthquake characterization using single-station waveforms and a convolutional neural network. *Seismol. Res. Lett.* **90**, 517–529 (2019).

37. Zhang, X., Zhang, M. & Tian, X. Real-time earthquake early warning with deep learning: application to the 2016 Central Apennines, Italy earthquake sequence. *Geophys. Res. Lett.* **48**, 2020GL089394 (2021).

38. Hayes, G. P. et al. Slab2, a comprehensive subduction zone geometry model. *Science* **362**, 58–61 (2018).

39. LeCun, Y., Bottou, L., Bengio, Y. & Haffner, P. Gradient-based learning applied to document recognition. *Proc. IEEE* **86**, 2278–2324 (1998).

40. Bletyer, Q. et al. A detailed source model for the M<sub>w</sub> 9.0 Tohoku-Oki earthquake reconciling geodynamics, seismology, and tsunami records. *J. Geophys. Res. Solid Earth* **119**, 7636–7653 (2014).

41. Bóse, M. et al. FinDer v2: improved real-time ground-motion predictions for M2–M9 with seismic finite-source characterization. *Geophys. J. Int.* **212**, 725–742 (2018).
42. Fukushima, K. & Miyake, S. Neocognitron: a new algorithm for pattern recognition tolerant of deformations and shifts in position. *Pattern Recognit.* **15**, 455–469 (1982).

43. Rumelhart, D., Hinton, G. & Williams, R. Learning representations by back-propagating errors. *Nature* **323**, 533–536 (1986).

44. Li, Z., Zhang, Y. & Arora, S. Why are convolutional nets more sample-efficient than fully connected nets? Preprint at https://arxiv.org/abs/2010.08515 (2021).

45. Krizhevsky, A., Sutskever, I. & Hinton, G. E. ImageNet classification with deep convolutional neural networks. In *Advances in Neural Information Processing Systems 25 (NIPS 2012)* (eds Pereira, F. et al.) (2012).

46. Pham, H., Dai, Z., Xie, Q. & Le, Q. V. Meta pseudo labels. In *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition* 11557–11568 (2021).

47. Bergen, K. J., Johnson, P. A., de Hoop, M. V. & Beroza, G. C. Machine learning for data-driven discovery in solid Earth geoscience. *Science* **363**, eaau0323 (2019).

48. Ren, C. X., Hulbert, C., Johnson, P. A. & Rouet-Leduc, B. Machine learning and fault rupture: a review. In *Machine Learning in Geosciences* (eds Moseley, B. & Krischer, L.) 57–107 (Elsevier, 2020).

49. Kingma, D. P. & Ba, J. L. Adam: a method for stochastic optimization. In *Intl Conf. Learning Representations* (eds Bengio, Y. & LeCun, Y.) (ICLR, 2015).

50. Zhu, W., Mousavi, S. M. & Beroza, G. C. Seismic signal denoising and decomposition using deep neural networks. *IEEE Trans. Geosci. Remote Sens.* **57**, 9476–9488 (2019).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_60.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022
# Methods

## Construction of the training database

We train PEGSNet on a database of synthetic data augmented with empirical noise. For each example in the database, three-component waveforms at each station location are obtained as follows.

**Source parameters.** The event source location is randomly picked from 1,400 possible locations extracted from the Slab2.0 model of subduction geometry³⁸ at two different depths (20 and 30 km). Given latitude, longitude and depth, strike and dip angles are determined by the subduction geometry and rake is randomly extracted from a normal distribution with mean = 90° and standard deviation of 10°. The event final $M_w$ is drawn from a uniform distribution with min = 5.5 and max = 10.0. We deliberately choose not to use a Gutenberg-Richter distribution for $M_w$ to avoid a sampling bias during training, by which the model might better estimate certain magnitude values simply because they are more represented in the training database. Finally, from $M_w$ we compute the scalar moment $M_0 = 10^{1.5M_w+9.1}$.

**Source time function.** Given $M_0$, a pseudo-empirical STF is computed using the STF model described in a previous work⁵, which includes a multiplicative error term and is valid for earthquakes with $M_w > 7.0$. In summary,

$$STF(t) = M_0 \frac{f(t)}{\int f(t) dt}, \qquad (1)$$

with:

$$f(t) = t \exp\{-0.5(\lambda t)^2\} [1 + N(t)], \qquad (2a)$$

$$\lambda = 10^{(7.24-0.41\log(M_0)+\epsilon)}, \qquad (2b)$$

$$N(t) = 0.38 \frac{n(t)}{\sigma}, \qquad (2c)$$

where $\epsilon$ is drawn from a Gaussian distribution with zero mean and standard deviation of 0.15, $n(t)$ is the time integral of a Gaussian noise time series with zero mean and $\sigma$ is the standard deviation of $n(t)$. The term $\epsilon$ accounts for variability in the STF duration for a given $M_0$, whereas $N(t)$ models the characteristics of noise observed in real STFs⁵. Examples of final STFs for different magnitude values are shown in Extended Data Fig. 1.

**Computing synthetic waveforms.** With the selected source parameters and STFs, we use the normal-mode approach described in a previous work¹⁴ to compute three-component synthetic waveforms in a spatial domain of about 20° around the source epicentre. The resulting seismometer responses are convolved with the STF of the corresponding synthetic event and multiplied by the scalar moment to obtain synthetic traces of acceleration sampled at 1 Hz at each station location. Finally, traces are bandpass-filtered between 2.0 mHz (Butterworth, two poles, causal) and 30.0 mHz (Butterworth, six poles, causal). The final seismograms are 700 s long centred at the event origin time.

**Noise database.** The noise database consists of 259 days of three-component waveforms for two non-continuous time intervals: between January 2011 and October 2011 (excluding March 2011) and between January 2014 and April 2014. These intervals have been chosen to sample variable (seasonal) noise conditions. We note that the temporal range spanned by the noise database does not overlap with any of the earthquakes used for real data cases (Extended Data Fig. 10a). We first divide the daily recordings into 1-h-long traces and then apply limited preprocessing, removing the instrument response, the mean and the

linear trend, converting to acceleration and decimating the original traces from 20 to 1 Hz. Finally, each trace is filtered using the same bandpass filter applied to the synthetic seismograms (see previous step) and stored. Note that no a priori assumptions on levels and characteristics of the selected noise are made. On the contrary, we include all real noise conditions found in continuous seismic recordings in the specified period range. This is because, in principle, we want the model to be able to generalize well under a broad range of noise conditions.

**Adding empirical noise to synthetics.** From the noise database described in the previous step, a realization of noise (700 s long) is extracted by randomly selecting a starting time point. In this process, we make sure to use different noise data in training, validation and test sets. To preserve spatial coherence of noise across the seismic network, the same time interval is used for all stations for a given event. The selected noise traces are then added to the corresponding acceleration seismograms to produce the final input data for PEGSNet. If noise data are not available for one or more stations in the selected time window for a given event, we discard those stations by setting the corresponding final trace amplitudes (noise and PEGS) to zero in the input data.

**Preprocessing of input data.** Before being fed to PEGSNet, we first sort the input waveform for each example based on station longitude. We found this approach to be effective, but we note that the problem of concatenating station waveforms in a meaningful way in deep learning is an active area of research³⁵. Then, on the basis of the theoretical P-wave arrival time ($T_p$) at each station for a given event, we set the amplitude of the seismograms to zero for $t \ge T_p$. Note that PEGSNet does not perform P-wave triggering itself. Instead, it relies on theoretical P-wave arrivals. In a real-time scenario, any existing P-wave triggering algorithm (whether based on machine learning or not) can be used to set the switch on the corresponding stations whose data can then be passed to PEGSNet.

To limit the influence of very noisy traces and to suppress high amplitudes (possibly related to the background regional seismicity), we further clipped the resulting trace using a threshold of ±10 nm s⁻². This threshold is chosen according to the maximum PEGS amplitude for an $M_w = 10$ earthquake at 315 s as found in the calculated database of noise-free synthetics. Amplitudes are finally scaled by 10 nm s⁻² to facilitate convergence of the optimizer, and at the same time, to preserve information about the relative amplitudes of the PEGS radiation pattern across the seismic network. Finally, to simulate missing data and/or problematic sensors, we randomly mute 5% of the stations for each event by setting to zero the amplitudes of the corresponding traces.

## Deep learning and PEGSNet

**Previous work.** Convolutional neural networks (CNNs) originated in neocognitron⁴² and became practical once it was found that the backpropagation procedure⁴³ can be used to compute the gradient of an objective function with respect to the weights of the network. CNNs are a regularized form of neural networks, that is, the function space they represent is simpler and they are more sample-efficient than fully connected neural networks⁴⁴. Deep CNNs have brought about a revolution in computer vision, and have had a role in almost every state-of-the-art approach for tasks related to recognition and detection in images⁴⁵,⁴⁶. In geoscience, machine learning has shown strong potential for data-driven discovery of previously unknown signals and physical processes hidden in large volumes of noisy data⁴⁷,⁴⁸.

We note, however, that our choice of a deep learning model over classical machine learning models offers an appealing framework to directly deal with raw seismograms. As a consequence, this choice enables us to explore a larger function space that is not limited by building an a priori set of features, which is a requirement for applying classical machine learning models on seismogram data.
Article

Successful applications of deep learning in seismology have provided new tools for pushing the detection limit of small seismic signals<sup>31,32</sup> and for the characterization of earthquake source parameters (magnitude and location)<sup>33-35</sup> with EEWS applications<sup>29,36,37</sup>. We present a deep learning model, PEGSNet, trained to estimate earthquake location and track the time-dependent magnitude, *M*<sub>w</sub>(*t*), from PEGS data, before P-wave arrivals.

**Description of PEGSNet architecture.** PEGSNet is a deep CNN that combines convolutional layers and fully connected layers in sequence (Extended Data Fig. 2a). The input of the network is a multi-channel image of size (*M*, *N*, *c*) where *M* is 315 (corresponding to 315-s-long traces sampled at 1 Hz), *N* is the number of stations (74) and *c* is the number of seismogram components used (three: east, north and vertical). The outputs of the network are three values corresponding to moment magnitude (*M*<sub>w</sub>), latitude (φ) and longitude (λ), where *M*<sub>w</sub> is time dependent. The training strategy used to learn *M*<sub>w</sub>(*t*) from the data is described below.

The first part of the model (the CNN) consists of eight convolutional blocks. Each block is made of one convolutional layer with a rectified linear unit (ReLU) activation function followed by a dropout layer. The number of filters in each convolutional layer increases from 32 (blocks 1–5) to 64 (blocks 6–7) to 128 (block 8) to progressively extract more detailed features of the input data. A fixed kernel size of 3 × 3 is used in each convolutional layer. We use spatial dropout with a fixed rate of 4% to reduce overfitting of the training set. Maximum pooling layers are added starting from block 4 to reduce the overall dimension of the input features by a factor of 4. The output of the CNN is then flattened and fed to a sequence of two dense layers of size 512, and 256 with a ReLU activation function and standard dropout with a 4% rate. Fully connected layers perform the high-level reasoning and map the learned features to the desired outputs. The output layer consists of three neurons that perform regression through a hyperbolic tangent activation function (tanh). The labelling strategy for *M*<sub>w</sub>(*t*), φ and λ is discussed in detail below. The total number of parameters in the network is 1,479,427.

**Learning strategy.** The purpose of the model is to track the moment released by a given earthquake as it evolves from the origin time. A specific learning strategy has been developed to address this task (Extended Data Fig. 2).

**Labelling.** Labels are φ, λ and a time-dependent *M*<sub>w</sub>, φ, λ simply correspond to the true values for each event. *M*<sub>w</sub>(*t*) is the time integration of the STF for each event. As detailed in the next section, the model is trained by randomly perturbing the ending time of the input seismograms, so that for a given ending time the input data are associated with the value of *M*<sub>w</sub>(*t*) at that time. To enforce the role of the tanh activation function in the output layer, we further scale all the labels to fall in the [-1, 1] interval through min/max normalization.

**Learning the time-dependent moment release.** In order for PEGSNet to learn *M*<sub>w</sub>(*t*), we randomly perturb the starting time of the input data during training (Extended Data Fig. 2c). Every time that an example is extracted from the dataset, a value (*T*<sub>1</sub>) is drawn at random from a uniform distribution between -315 and 0 (s). *T*<sub>1</sub> is the time relative to the earthquake origin time (*T*<sub>0</sub>) corresponding to the starting time of the selected seismograms for that example. In practice, from the 700-s-long seismograms (centred on *T*<sub>0</sub>) in the database, we extract traces from *T*<sub>1</sub> to *T*<sub>2</sub> = *T*<sub>1</sub> + 315 s: for *T*<sub>1</sub> = -315 s the extracted traces end at *T*<sub>0</sub>; for *T*<sub>1</sub> = 0 s the traces start at *T*<sub>0</sub> and end 315 s after. Once a value for *T*<sub>1</sub> is selected, the value of *M*<sub>w</sub>(*T*<sub>2</sub>) is assigned as the corresponding label for this example (Extended Data Fig. 2d). This enables the model to learn patterns in the data as the STF evolves with time.

**Training.** The full database (500,000 examples of synthetic earthquakes) is split into training (350,000) validation (100,000) and test (50,000) sets, following a 70/20/10 strategy. The network is trained for

200 epochs (using batches of size 512) on the training set by minimizing the Huber loss between the true values and the predicted earthquake source parameters using the Adam algorithm<sup>49</sup>, with its default parameters (β<sub>1</sub> = 0.9 and β<sub>2</sub> = 0.999) and a learning rate of 0.001. At the end of each epoch, the model is tested on the validation set to assess the learning performance and avoid overfitting (Extended Data Fig. 2b). After learning, the model that achieved the best performance (lowest loss value) on the validation set is selected as the final model. The final model is then tested against the test set (therefore with data that has never been seen by PEGSNet during training) to assess its final performance.

**Testing strategy.** Once PEGSNet is trained, it can be used to estimate *M*<sub>w</sub>(*t*) in a real-time scenario. We assess the latency performance of PEGSNet on the test set with the following procedure (Extended Data Fig. 4). For each example in the test set, we slide a 315-s-long window [*T*<sub>1</sub>, *T*<sub>2</sub> = *T*<sub>1</sub> + 315 s] through the data with a time step of 1 s. The starting window ends at the earthquake origin time *T*<sub>0</sub> (*T*<sub>2</sub> = *T*<sub>0</sub> and *T*<sub>1</sub> = *T*<sub>0</sub> - 315 s) and the final window starts at the earthquake origin time (*T*<sub>2</sub> = *T*<sub>0</sub> + 315 s and *T*<sub>1</sub> = *T*<sub>0</sub>). We let PEGSNet predict *M*<sub>w</sub>(*T*<sub>2</sub>) at each time step, thus progressively reconstructing the STF. Each *M*<sub>w</sub>(*t*) estimate made by PEGSNet only makes use of information in the past 315 s. The same procedure is also applied to real data (Fig. 3 and Extended Data Fig. 10), to simulate a playback of the data as if they were fed to PEGSNet in real-time.

**Test on noise-free synthetics.** We investigate the performance of PEGSNet by using the same database described above but without including noise in the input data. Training and testing on noise-free synthetic data provides an upper limit of PEGSNet performance. Although this experiment represents a virtually impossible scenario for real-world applications, the results can reveal inherent limitations of our model or in the input data. Extended Data Fig. 6a shows the accuracy map for the test set. As expected, the model is able to determine the final *M*<sub>w</sub> of the events with high accuracy and similar performance regardless of the actual *M*<sub>w</sub> of the event, except at early times. To look at the latency performance more in detail, Extended Data Fig. 6b shows the density plot of the residuals as a function of time for the whole noise-free test set. Errors are mostly confined within ±0.1 but are relatively higher in the first 10–15 s after origin. We relate this observation to a combination of two factors: first, in the first 15 s after origin, very small PEGS amplitudes are expected at very few stations in the seismic network, partially owing to the cancelling effect between the direct and induced terms. This can lead to a situation in which little information is present in the input images and the model ends up predicting the mean value of the labels at these early times. Second, the seismic network geometry may not be optimal for recording PEGS amplitudes in this time window. Finally, we note that similar behaviour is observed for the results obtained on the noisy database (Fig. 2a, b) but with a higher latency (30–40 s). This highlights the role of noise in degrading the optimal performance of PEGSNet.

Preprocessing of real data

The preprocessing steps for real data closely follow the procedure detailed in previous work<sup>7</sup>. For each station and each component:

1. Select 1-h-long raw seismograms ending at the theoretical *T*<sub>p</sub> calculated using the source location from the United States Geological Survey (USGS) catalogue (Extended Data Fig. 10a);

2. Remove the mean;

3. Remove the instrument response and obtain acceleration signals;

4. Lowpass 30.0 mHz (Butterworth, six poles, causal);

5. Decimate to 1 Hz;

6. Highpass 2.0 mHz (Butterworth, two poles, causal);

7. Clip to ±10 nm s⁻² and scale by the same value;

8. Pad with zeros for *t* ≥ *T*<sub>p</sub> and select 700-s-long trace centred on *T*<sub>0</sub>.
This procedure is the same as that used to generate the synthetic database, except that here, traces need to be cut at the P-wave arrival first to avoid contamination of PEGS by the P-wave during instrument response removal.

To speed up our testing procedure (see Methods subsection ‘Testing strategy’), the data are preprocessed once and then sliced into input for PEGSNet at each time step. In an online version of our model, this is unfeasible as all the preprocessing steps need to be applied each time that new packets of data are streamed in. We simulate the conditions of a real-time workflow on the Tohoku-Oki data to assess potential discrepancies with the simplified workflow in the results: at each time step, we apply the preprocessing steps described above, using 1-h-long trace ending at the current time step. We find that the resulting PEGSNet predictions obtained using the two workflows are essentially indistinguishable from each other (Extended Data Fig. 11).

## Predictions on additional real data

We tested PEGSNet on all the subduction earthquakes (dip-slip mechanism within 40 km from the megathrust) with $M_w \ge 7$ that have occurred since January 2003, without considering aftershocks (Extended Data Fig. 10). Among them, the 2003 $M_w = 8.2$ Hokkaido earthquake is at the edge of PEGSNet’s lower sensitivity limit of 8.3. For this event, PEGSNet estimates the final $M_w$ after about two minutes (Extended Data Fig. 10b), in agreement with what was previously observed on the test set for events with similar magnitude (Fig. 2a). However, given the expected lower accuracy and higher errors for this event, we consider these predictions less reliable. For lower-magnitude events, PEGSNet predictions converge toward the noise baseline of 6.5 or never exceed its lower sensitivity limit, confirming that PEGS from $M_w < 8.0$ earthquakes are essentially indistinguishable from noise (Extended Data Fig. 10c–f). Deep learning denoising techniques for coherent noise removal⁵⁰ might prove successful in enhancing PEGSNet performance and will be the subject of future work.

## Data availability

Waveforms and metadata used in this study were provided by the Global Seismograph Network⁵¹, the New China Digital Seismograph Network⁵², GEOSCOPE⁵³ (IPGP/EOST, 1982), the Japan Meteorological Agency Seismic Network (https://www.jma.go.jp/en/quake) and the National Research Institute for Earth Science and Disaster Resilience/F-net⁵⁴, all publicly available from the IRIS Data Management Center (https://ds.iris.edu/ds/nodes/dmc) and the F-net data centre (https://www.fnet.bosai.go.jp/top.php). Source data are provided with this paper.

## Code availability

The codes used to generate and analyse the data shown within this paper are available from the corresponding author upon request. PEGSNet is built and trained using PyTorch⁵⁵. Waveforms were retrieved from IRIS with Obspy⁵⁶ and obspyDMT⁵⁷. Figures were produced with the Generic Mapping Tool (GMT)⁵⁸ and matplotlib⁵⁹.

51. Albuquerque Seismological Laboratory (ASL)/USGS. Global Seismographic Network (GSN-IRIS/USGS) (International Federation of Digital Seismograph Networks, 1988); https://doi.org/10.7914/SN/IU.

52. ASL/USGS. New China Digital Seismograph Network (NCDSN) (International Federation of Digital Seismograph Networks, 1992); https://doi.org/10.7914/SN/IC.

53. Institut de physique du globe de Paris (IPGP) and Ecole et Observatoire des Sciences de la Terre de Strasbourg (EOST). GEOSCOPE, French Global Network of Broad Band Seismic Stations (IPGP, 1982); https://doi.org/10.18715/GEOSCOPE.G.

54. Okada, Y. et al. Recent progress of seismic observation networks in Japan — Hi-net, F-net, K-NET and KiK-net. Earth Planets Space 56, xv–xxviii (2004).

55. Paszke et al. PyTorch: an imperative style, high-performance deep learning library. In Advances in Neural Information Processing Systems 32 (NeurIPS 2019) (eds Wallach, H. et al.) 8024–8035 (Curran Associates, 2019).

56. Beyreuther, M. et al. Obspy: a Python toolbox for seismology. Seismol. Res. Lett. 81, 530–533 (2010).

57. Hosseini, K. & Sigloch, K. ObspyDMT: a Python toolbox for retrieving and processing large seismological data sets. Solid Earth 8, 1047–1070 (2017).

58. Wessel, P. et al. The Generic Mapping Tools version 6. Geochem. Geophys. Geosyst. 20, 5556–5564 (2019).

59. Hunter, J. D. Matplotlib: a 2D graphics environment. Comput. Sci. Eng. 9, 90–95 (2007).

**Acknowledgements** This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement 949221). This work has been supported by the French government, through the UCAJEDI Investments in the Future project managed by the National Research Agency (ANR), ANR-15-IDEX-01. This work was granted access to the HPC resources of IDRIS under the allocations 2020-AD011012142, 2021-AP011012536 and 2021-A0101012314 made by GENCI. B.R.-L.’s work was supported by Institutional Support (LDRD) at Los Alamos (20200278ER). Numerical computations for the synthetic database of PEGS were performed on the S-CAPAD platform, IPGP, France. We thank M. Böse for providing the results of FinDer2 on Tohoku-Oki data.

**Author contributions** A.L. developed PEGSNet and produced all the results and figures presented here. Q.B. had the original idea and supervised the study. B.R.-L. designed the preliminary algorithm and provided the raw noise data. J.-P.A. provided expertise in modelling PEGS and realistic source time functions. K.J. provided the synthetic database of PEGS necessary to produce the training database. All the authors contributed to writing the manuscript.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41586-022-04672-7.

**Correspondence and requests for materials** should be addressed to Andrea Licciardi.

**Peer review information** Nature thanks Benjamin Fong Chao and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
<figure><img src="image_61.png" /></figure>

**Extended Data Fig. 1** | STF database. a, STFs (moment rate) coloured by their final magnitude in the training database from a previous model⁵ (see description in Methods). Only STFs for events with final *M*<sub>w</sub> between 7.0 and 9.5 are shown. b, The complete *M*<sub>w</sub>(*t*) distribution is used for labelling the database.
<figure><img src="image_62.png" /></figure>

**Extended Data Fig. 2** | PEGSNet architecture and training strategy.
a, The input data for one example consists of a three-channel image of shape $M \times N$, where $M$ is the number of time samples and $N$ is the number of seismic stations. Only the vertical component of the input data is displayed for simplicity. Each convolutional block is composed of a convolutional layer (yellow) with a ReLu activation (orange), a spatial dropout layer (light blue). Max pooling layers (red) reduce each dimension of the input data by a factor of two. The number of filters used in each convolutional layer is indicated for clarity. The last convolutional block is connected to dense layers (purple) (using a ReLu activation function), with dropout (light blue). The output layer uses a tanh activation function to predict values of $M_w$, latitude ($\phi$) and longitude ($\lambda$).
b, The value of the Huber loss is plotted as a function of epochs

for the training (Train.) and validation (Val.) sets. Each epoch corresponds to a full pass over the training set in batches of size 512. The red star indicates the epoch with the minimum value of the loss on the validation set. The corresponding model is used for predictions on the test set and real data.
c, Data from one example from the training database (vertical component). The grey shaded area corresponds to the input data for PEGSNet shown in a. $T_1$ and $T_2$ are the beginning and end of the selected input window. During training, $T_1$ is selected at random and $T_2 = T_1 + 315$ s. d, Moment rate (blue) and $M_w(t)$ (dark grey) for the selected event. Given the randomly selected value for $T_1$ for this example, the corresponding label is $M_w(T_2)$, that is, at the end of the selected window. This is compared with the predicted $M_w$ estimated by PEGSNet in a and used for training.
<figure><img src="image_63.png" /></figure>

<figure><img src="image_64.png" /></figure>

<figure><img src="image_65.png" /></figure>

<figure><img src="image_66.png" /></figure>

<figure><img src="image_67.png" /></figure>

<figure><img src="image_68.png" /></figure>

**Extended Data Fig. 3** | Results of predictions on the test set for location. **a, b,** Analysis of the residuals as a function of time (only the first 180 s from origin time are shown) for latitude (**a**) and longitude (**b**). Boxes correspond to the third to first (Q3–Q1) interquartile range, the black line within the boxes is the median and whiskers indicate the 5th and 95th percentiles. **c, d,** Density

plot of true versus predicted latitude (**c**) and longitude (**d**) values at 180 s after origin (red boxes in **a, b**). **e, f,** Corresponding histogram of the residuals with reported mean and standard deviations for latitude (**e**) and longitude (**f**) at 180 s after origin.
<figure><img src="image_69.png" /></figure>

**Extended Data Fig. 4 | Testing procedure.** Each input example (top row) is parsed with a time window of length 315 s sliding with 1-s time step. At each time step, the data between $T_1$ and $T_2$ are passed to PEGSNet (middle row). For display purposes, only four specific time steps are shown here with their $T_1$ and $T_2$ indicated for clarity. The window with $T_1$ = -200 s and $T_2$ = 115 s corresponds to

the grey shaded area in the top row. At each time step, PEGSNet makes a prediction of $M_w$ at the end of the input window, that is, $M_w(T_2)$ (blue circles in the bottom figure), to reconstruct the STF (yellow line) in a real-time scenario. Red dashed lines indicate the value of $T_2$ of each input window.
<figure><img src="image_70.png" /></figure>

**Extended Data Fig. 5** | Effect of noise on PEGSNet predictions. a, Frequency distribution of the mean standard deviation of the noise ($\sigma$) for the whole seismic network, in a pre-event time window of five minutes, for test set only. The solid red line is the median of the distribution (second quartile, Q2), dashed lines are the first (Q1) and third (Q3) quartiles. b–e, As in Fig. 2a, but computed

on subsets of the test set for which $\sigma < Q1$ (b), $Q1 < \sigma < Q2$ (c), $Q2 < \sigma < Q3$ (d), and $\sigma > Q3$ (e). Note the high accuracy below 8.3 after 150 s from origin time for low-noise conditions (b). These maps give empirical limits on the noise levels that could potentially be useful to interpret the performance of our model in real time.
<figure><img src="image_71.png" /></figure>

**Extended Data Fig. 6** | Results on noise-free test set. a, As in Fig. 2a, but for a test set including only noise-free PEGS waveforms. b, Density plot of the residuals as a function of time for the same test set. The solid red line is the

median, the dashed lines are the 5th and 95th percentiles of the distribution.
Note that both in a and b, the predictions are obtained for a model that has been
trained on a database of noise-free synthetics.
<figure><img src="image_72.png" /></figure>

**Extended Data Fig. 7** | Comparison of PEGSNet with existing EEWS. a, b, As in Fig. 3, but with two additional EEWS algorithms included. The results for Tohoku of finite-fault GNSS-based G-larmS²⁸ and seismic point-source ElarmS

are taken from an earlier work²⁸. In the main text, only the best-performing algorithm for each type (JMA, point-source; FinDer2, seismic finite-fault; and BEFORES, GNSS-based finite-fault) are shown.
<figure><img src="image_73.png" /></figure>

**Extended Data Fig. 8 | Synthetic test with Tohoku STF.** Density plot of the predictions on 1,000 examples obtained by combining different noise recordings (randomly extracted from the test set) with PEGS waveforms, from a synthetic source that mimics that of the Tohoku-Oki earthquake. The median, Q1-Q3 interquartile range and the 5th-95th percentiles of the distribution are reported with solid, dotted and dashed red lines, respectively. The synthetic Tohoku-like source is obtained from the 'true' STF (orange line, same as in Fig. 3 with ±0.3 *M*<sub>w</sub> indicated by dashed lines), the hypocentre location (Extended Data Fig. 10a), and the following values for the strike, dip and rake: 193.0°, 8.9° and 78.4° (USGS catalogue). The workflow described in Methods is used to

combine synthetics and noise. PEGSNet prediction for real Tohoku data is reported for comparison (blue line, same as in Fig. 3). Even if it has been trained on synthetic data plus noise, PEGSNet is able to generalize well to real data (the blue line is within the expected variability of the predictions obtained on synthetics). In the first 50 s after origin, the variability of the predictions is high and strongly affected by noise (similar to Fig. 2c), whereas after 50 s, predictions have similar *M*<sub>w</sub>(t) values and therefore are more robust. Note that for the lowest noise conditions, faster PEGSNet response is virtually possible, as early as about 30 s after origin (as indicated by the 95th percentile of the distribution).
<figure><img src="image_74.png" /></figure>

**Extended Data Fig. 9 | PEGSNet noise test.** Density plot of the predictions obtained for 1,000 real noise recordings (the same used in Extended Data Fig. 8) assuming P-wave arrivals (but no PEGS) as in the real data for the Tohoku-Oki earthquake. The median, Q1–Q3 interquartile range and the 5th–95th percentiles of the distribution are reported with solid, dotted and dashed

red lines respectively. This test confirms that the results of Fig. 3 are indeed constrained by the data in the pre-P-wave time window (PEGS). The predicted $M_w(t)$ tends towards a constant value of 6.5, which is considered as the baseline value when no information can be extracted from PEGS.
a

<table><thead><tr><th>Name</th><th>M<sub>w</sub></th><th>Longitude</th><th>Latitude</th><th>Depth [km]</th></tr></thead><tbody><tr><td>2003-09-25T19:50:06</td><td>8.2</td><td>143.910</td><td>41.815</td><td>27.0</td></tr><tr><td>2003-10-31T01:06:28</td><td>7.0</td><td>142.619</td><td>37.812</td><td>10.0</td></tr><tr><td>2005-08-16T02:46:28</td><td>7.2</td><td>142.039</td><td>38.276</td><td>36.0</td></tr><tr><td>2008-07-19T02:39:28</td><td>7.0</td><td>142.214</td><td>37.552</td><td>22.0</td></tr><tr><td>2011-03-09T02:45:20</td><td>7.3</td><td>142.842</td><td>38.435</td><td>32.0</td></tr><tr><td>2011-03-11T05:46:24</td><td>9.1</td><td>142.373</td><td>38.297</td><td>29.0</td></tr></tbody></table>

b

<figure><img src="image_75.png" /></figure>

c

2003/10/31 01:06:28 -- Off East Coast of Honshu

<figure><img src="image_76.png" /></figure>

d

2005/08/16 02:46:28 -- Near East Coast of Honshu

<figure><img src="image_77.png" /></figure>

e

2008/07/19 02:39:28 -- Off East Coast of Honshu

<figure><img src="image_78.png" /></figure>

f

2011/03/09 02:45:20 -- Near East Coast of Honshu

<figure><img src="image_79.png" /></figure>

**Extended Data Fig. 10 | PEGSNet predictions on real data.** a, List of earthquakes used for the real data cases. All parameters are taken from the USGS catalogue. Earthquakes were selected according to the following criteria: $M_w \ge 7.0$, dip-slip focal mechanism and located within 40 km from the source locations used for training (green lines in Fig. 1). **b–f**, PEGSNet $M_w(t)$ predictions are indicated by the blue line. The integrated STFs (orange lines) are taken from the SCARDEC database with dashed lines representing $\pm 0.3$ magnitude units.

The red horizontal line marks the empirical lower limit of PEGSNet sensitivity (8.3). Note that for the 2003 $M_w$ 8.2 Hokkaido earthquake (b), this limit is crossed several times after 120 s from origin time. Out of four events with final $M_w$ below 8 (c–f), three (d–f) show a predicted $M_w(t)$ that is constant around the noise baseline (see Extended Data Fig. 9) indicating that no information can be extracted from the pre-P-wave time window.
<figure><img src="image_80.png" /></figure>

**Extended Data Fig. 11 | PEGSNet predictions for two data workflows.**

The predictions for the Tohoku-Oki data obtained through the original workflow described in Methods (blue line, same as in Fig. 3) are compared with those obtained with a modified workflow (red line). Instead of preprocessing the data once and then slice them into input windows for PEGSNet (original

workflow), the data are preprocessed and fed to PEGSNet at each time step (modified workflow, see details in Methods) to simulate an online scenario. The results of two workflows show negligible differences. Note that the blue line is twice as thick as the red line. For reference, the orange lines indicate the true STF (solid) and ±0.3 magnitude units (dashed).
## 3 Perspectives

Le développement de PEGSnet est un projet d'envergure, financé par l'ERC Starting EARLI. Ce financement me permet d'encadrer une équipe de 2 doctorants et 2 postdoctorants, Gabriela Arias, Pablo Lara, Andrea Licciardi et Kévin Juhel. À court terme, nous travaillons sur des applications régionales ainsi que sur une technique de débruitage des PEGS par IA. À moyen terme, nous envisageons une application globale : un algorithme pouvant estimer quasi-instantanément la magnitude et la localisation des grands séismes partout sur Terre à partir de données sismologiques disponibles en temps réel partout sur Terre. Enfin, nous travaillons sur un algorithme d'IA complémentaire pour l'alerte des séismes de magnitude inférieure à 8.

### 3.1 Applications régionales

#### Rétrospective du séisme de Maule (M<sub>w</sub> 8.8, Chili 2010)

Dans le cadre de sa thèse débutée en Janvier 2021, Gabriela Arias exporte au Chili (Figure IV.2) l'approche développée par Andrea Licciardi au Japon.

<figure><img src="image_81.png" /><figcaption>**Figure IV.2** – Stations ayant enregistré le séisme de Maule et sources synthétiques considérées pour construire la base de données d'entrainement (équivalent de la Figure 1 de la partie IV.2 pour application au Chili) [Arias et al., in prep.].</figcaption></figure>
Bien que les données soient beaucoup plus bruitées et moins nombreuses qu'au Japon, les résultats indiquent que PEGSnet aurait pu estimer la magnitude du séisme de Maule après une centaine de secondes (Figure IV.3). La magnitude minimale pour que PEGSnet fonctionne dans la configuration du séisme de Maule est toutefois $M_w = 8.6$ – en dessous de cette valeur PEGSnet échoue à discerner les PEGS du bruit – [Arias et al., in prep.], très au dessus du seuil dans la configuration du séisme de Tohoku-Oki.

<figure><img src="image_82.png" /><figcaption>**Figure IV.3** – Estimation de la magnitude du séisme de Maule (courbe bleue) en fonction du temps comparée à l'évolution réelle de la magnitude (solution SCARDEC [Vallée et al., 2011]) et à la même estimation pour 1,000 échantillons de bruit empirique sans PEGS (courbes grises) [Arias et al., in prep.].</figcaption></figure>

## Implémentation dans le système d'alerte Péruvien

Les études évoquées ci-avant démontrent le potentiel de PEGSnet pour l'alerte rapide, en particulier tsunami. L'étape suivante est l'implémentation concrète dans un système d'alerte opérationnel. En collaboration avec Adolfo Inza à l'Instituto Geofísico del Peru (IGP), nous travaillons à l'implémentation de PEGSnet en temps réel sur le réseau de stations large-bande, riche d'une soixantaine de stations dans lequel le gouvernement péruvien a récemment investi (Figure IV.4). Cette application concrète au Pérou est une priorité à moyen terme, avec des visites longues durées dans les 2 sens. J'ai ainsi passé 2 mois à l'IGP en 2021, accompagné
d'Andrea Licciardi et Gabriela Arias. Adolfo Inza fera le chemin inverse pour 2 mois en 2022.

<figure><img src="image_83.png" /><figcaption>**Figure IV.4** – Stations large-bande transmettant en temps-réel au Pérou (triangles verts) et séismes importants ($M_w$ > 7.5) depuis 1950.</figcaption></figure>

## 3.2 Débruitage des PEGS par IA

Nous avons vu dans la section IV.2 (Extended Data Figure 6) que si les données étaient dépourvues de bruit, PEGSnet permettrait une estimation instantanée de la magnitude de n'importe quel séisme, quelque soit sa magnitude. PEGSnet n'est donc limité que par le bruit sur les stations. Si l'on parvenait à éliminer ce bruit, on pourrait donc améliorer très significativement les performances de PEGSnet, en particulier pour les séismes de magnitude modérée ($7 < M_w < 8.5$). Nous développons ainsi un algorithme de débruitage des PEGS par Deep Learning. L'architecture de l'algorithme est très semblable à PEGSnet. Les données d'entrée sont les mêmes, mais au lieu d'estimer, en sortie, la latitude, la longitude et la magnitude de l'évènement, cet algorithme de débruitage estime les PEGS débruités (Figure IV.5).

L'algorithme prend avantage de la cohérence spatio-temporelle entre les PEGS attendus d'un bout à l'autre du réseau pour reconstruire les signaux et éliminer le bruit. Le résultat station par station est ainsi très supérieur à ce que l'on pourrait attendre d'un simple filtre
<figure><img src="image_84.png" /><figcaption>**Figure IV.5** – Exemple de débruitage des PEGS sur le réseau de stations de l'étude de la section IV.2. En haut les données observées, composées des PEGS et du bruit. En bas le signal attendu dépourvu du bruit. Au milieu, les signaux débruités par l'algorithme [Licciardi et al., in prep.].</figcaption></figure>

dans la gamme de fréquence des PEGS (Figure IV.6).

L'étape suivante est de configurer l'algorithme pour permettre un débruitage seconde par seconde, puis d'utiliser les PEGS débruités comme donnée d'entrée dans PEGSnet.

## 3.3 Vers un algorithme d'alerte global

Les caractéristiques exceptionnelles des PEGS permettent d'imaginer quelque chose d'in-
envisageable avec les signaux sismiques « classiques » : un algorithme d'alerte utilisant toutes
les stations de la planète transmettant en temps réel. En effet, les PEGS se propageant à la
vitesse de la lumière, il est, en théorie, possible d'enregistrer des signaux utiles aux antipodes
de l'épicentre avant l'arrivée du tsunami – et peut-être même avant les premières secousses
– dans les régions proches de la source. Un algorithme global de ce type présente plusieurs
intérêts. Le premier est la quantité d'information que l'on peut potentiellement extraire d'un
réseau global pour contraindre l'estimation. Le second est qu'un tel système pourrait fonction-
ner pour n'importe quel séisme partout sur Terre, offrant ainsi une « couverture » d'alerte dans
des régions ne bénéficiant pas de système opérationnel. Pour mettre au point un tel algorithme
<figure><img src="image_85.png" /><figcaption>**Figure IV.6** – Résultat station par station du débruitage sur le réseau entier [Licciardi et al., in prep.].</figcaption></figure>

global, l'obstacle majeur est de calculer les PEGS synthétiques sur la planète entière pour un ensemble exhaustif de paramètres de source. Construire la base de données synthétiques nécessaire à l'entrainement d'une version globale de PEGSnet constitue un challenge important sur le plan numérique. Dans le cadre de son postdoctorat, Kévin Juhel a réussi à construire cette base de données en adaptant une astuce numérique proposée par Zhang et al. [2020] afin de réduire drastiquement les temps de calcul [Juhel et al., in prep.]. L'étape suivante sera de collecter les données de toutes les stations transmettant en temps réel dans le monde afin de pouvoir réaliser l'entrainement.

## 3.4 Un algorithme alternatif efficace dans les toutes premières secondes

Nous avons montré que PEGSnet, notre algorithme d'IA basé sur les PEGS, est capable d'estimer de manière fiable et quasi-instantanée la magnitude des très grands séismes. Cela fait d'ores-et-déjà de PEGSnet un algorithme très utile pour l'alerte tsunami. La question demeure,
en revanche sur le potentiel de PEGSnet pour l'alerte sismique en temps que telle – i.e., la pré-
diction des mouvements forts. En l'état actuel – et même si nous travaillons sur des pistes pour
baisser ce seuil (section IV.3.2) – PEGSnet ne permet pas d'estimer la magnitude d'évènements
de magnitude inférieure à 8.3 (au Japon). La plupart des séismes n'atteignent une telle ma-
gnitude qu'après quelques dizaines de secondes, ce qui rendrait un système d'alerte n'utilisant
que PEGSnet inopérant durant ce temps. Pour combler ce manque, nous développons, dans le
cadre de la thèse de Pablo Espinoza Lara, un algorithme d'IA complémentaire qui estime la
magnitude et la localisation d'un séisme dans des temps records après son origine.

<figure><img src="image_86.png" /><figcaption>**Figure IV.7** – Performance de l’algorithme 3 secondes après l’arrivée de l’onde P [Lara et al., in prep.]. a) Densité de probabilité des estimations de la magnitude sur le jeu de données test (en abscisse la magnitude « réelle », en ordonnée l’estimation après 3 s). b) Idem pour la distance épicentrale. c) Idem pour la profondeur. d) Idem pour le back-azimut.</figcaption></figure>

La première phase de l’algorithme est la détection. Un premier algorithme d’IA opère conti-
nuellement, en temps réel, sur toutes les stations indépendamment. Lorsque un séisme est dé-
tecté sur une station, un deuxième algorithme est déclenché pour l’estimation. Cet algorithme
est entrainé sur une centaine de milliers de formes d’ondes compilée à partir du catalogue glo-
bal STEAD [Mousavi et al., 2019] et des catalogues régionaux du Chili, du Japon et du Pérou. L'objectif est d'estimer la magnitude, la distance épicentrale, le back-azimut et la profondeur du séisme à partir d'une seule station. L'entrainement n'est pas réalisé directement à partir des données brutes comme dans le cas de PEGSnet. 162 attributs sont calculés à partir des formes d'ondes dans les domaines temporel, spectral et cepstral et nourrissent la base de données d'entrainement. De manière intéressante, le domaine cepstral fournit les informations les plus importantes dans la décision de l'algorithme.

<figure><img src="image_87.png" /><figcaption>**Figure IV.8** – Performance de l'algorithme 30 secondes après l'arrivée de l'onde P [Lara et al., in prep.]. Idem à la figure IV.7 avec 30 s d'enregistrement après l'arrivée de l'onde P.</figcaption></figure>

La première estimation est réalisée à partir de 3 secondes d'enregistrement, i.e. 3 secondes après l'arrivée de l'onde P sur la première station. Les performances après 3 secondes indiquent qu'une estimation grossière de la magnitude est possible très tôt après le début de l'évènement (Figure IV.7.a), d'autant que l'estimation est réalisée sur une seule station, entrainant un gain de temps significatif par rapport aux approches multi-stations pour lesquels il est nécessaire d'attendre que le signal atteignent 4 ou 5 stations. Le back-azimut est également relativement
bien estimé après 3 s (Figure IV.7.d), autorisant une estimation grossière de la localisation, même si la distance épicentrale est moins bien estimée (Figure IV.7.b), à l'instar de la profondeur (Figure IV.7.c).

Avec 30 secondes d'enregistrement, les estimations sont toutefois nettement améliorées (Figure IV.8). Le progrès le plus spectaculaire est celui de la distance épicentrale (Figure IV.8.b), très fiablement estimée à partir de 30 s d'enregistrement sur une seule station. La magnitude est également mieux estimée (Figure IV.8.a), de même que la profondeur – même si celle-ci reste assez mal contrainte – (Figure IV.8.c). Aucun progrès significatif n'est observé pour le back-azimut. Globalement, la fiabilité des estimations augmente avec le nombre de secondes d'enregistrement fournies à l'algorithme (Figure IV.9), mais celui-ci montre un potentiel intéressant pour l'alerte très rapide.

<figure><img src="image_88.png" /><figcaption>**Figure IV.9** – Erreur sur les estimations de la magnitude (a), de la distance épicentrale (b), de la profondeur (c) et du back-azimut (d) en fonction du temps d'enregistrement après l'arrivée de l'onde P [Lara et al., in prep.].</figcaption></figure>

L'algorithme est implémenté en version bêta dans le système d'alerte péruvien et a montré des performances très encourageantes sur les séismes enregistrés ces derniers mois. Un grand nombre de stations accélérométriques transmettant en temps réel est actuellement déployé sur la côte péruvienne (tous les 30 km) dans le cadre du programme Sistema de Alerta Sísmica Peruano (SASPe). L'algorithme développé par Pablo y est systématiquement intégré.

Cet algorithme peut par ailleurs être utilisé pour d'autres applications que l'alerte sismique. Il est en effet capable de quantifier la sismicité à partir d'une seule station, ce qui peut s'avérer très utiles dans des régions mal instrumentées. Nous avons ainsi utilisé cet algorithme pour caractériser l'activité des répliques suivant le séisme de Nippes en Haïti (2021, M<sub>w</sub> 7.2) à partir d'une seule station à bas coût (très bruitée) installée chez l'habitant, avec des résultats équivalents, voire supérieurs, à ceux obtenus par les approches classiques sur le réseau régional
entier [Calais et al., 2022]. Des applications similaires sont planifiées, notamment sur l'île de
La Plata, en Équateur.
Conclusion

Nous n’avons pas répondu aux questions « Où ? » et « Quand ? » se dérouleront les prochains
grands séismes. Ces questions resteront probablement sans réponses pour de longues années
encore. Nous avons toutefois proposer des pistes pour avancer dans cette direction et apporté
des éléments pour mieux anticiper les séismes. Pourra-t-on un jour aller plus loin ?

La découverte récente des PEGS illustre qu’un signal sismologique majeur – le plus précoce
identifié à ce jour – est resté inconnu, caché dans le bruit sismique, pendant des années. L’idée
que d’autres signaux, plus précoces encore – i.e. précédant l’origine du séisme – restent en-
core à détecter dans le bruit sismique n’est donc pas déraisonnable. L’IA offre des perspectives
nouvelles pour extraire ce type de signaux. Pourra-t-elle, un jour, permettre la prédiction des
séismes ? La réponse dépend de la nature déterministe des ces phénomènes. Si les séismes sont
des processus purement chaotiques, se produisant sans phase précurserice, alors ils resteront
ontologiquement imprévisibles pour toujours. Si, en revanche, notre incapacité à prédire ces évè-
nements résulte de notre connaissance limitée des conditions physiques contrôlant la genèse des
ruptures en profondeur, alors la prédiction des séismes devient théoriquement possible. Même
si les conditions physiques en profondeur restent insuffisamment contraintes pour construire
des modèles capables de prédire de manière fiable les séismes à venir, l’IA est un outil pour
la détection « donnée orientée » qui ne requiert aucune modélisation ni aucune connaissance
des conditions physiques en profondeur. Si un signal précurseur existe dans les données, l’IA
représente ainsi un réel espoir de le découvrir.
Références bibliographiques

Ader, T. J., N. Lapusta, J.-P. Avouac, and J.-P. Ampuero (2014), Response of rate-and-state seismogenic faults to harmonic shear-stress perturbations, *Geophysical Journal International*, **198**(1), 385-413. 15

Aki, K., and P. G. Richards (2002), *Quantitative seismology*. 6

Allen, R. M. (2006), Probabilistic warning times for earthquake ground shaking in the san francisco bay area, *Seismological Research Letters*, **77**(3), 371-376. 84

Allen, R. M., and D. Melgar (2019), Earthquake early warning: Advances, scientific challenges, and societal needs, *Annual Review of Earth and Planetary Sciences*, **47**, 361–388. 21

Allgeyer, S., and P. Cummins (2014), Numerical tsunami simulation including elastic loading and seawater density stratification, *Geophysical Research Letters*, **41**(7), 2368-2375. 7

Altamimi, Z., P. Sillard, and C. Boucher (2002), Itrf2000: A new release of the international terrestrial reference frame for earth science applications, *Journal of Geophysical Research: Solid Earth*, 107(B10), ETG-2. 4

Altamimi, Z., X. Collilieux, J. Legrand, B. Garayt, and C. Boucher (2007), Itrf2005: A new release of the international terrestrial reference frame based on time series of station positions and earth orientation parameters, *Journal of Geophysical Research: Solid Earth*, 112(B9). 4

Altamimi, Z., X. Collilieux, and L. Métivier (2011), Itrf2008: an improved solution of the international terrestrial reference frame, *Journal of Geodesy*, 85(8), 457-473. 4
Altamimi, Z., P. Rebischung, L. Métivier, and X. Collilieux (2016), Itrf2014 : A new release of the international terrestrial reference frame modeling nonlinear station motions, *Journal of Geophysical Research : Solid Earth*, **121**(8), 6109–6131. 4

Amitrano, D. (2003), Brittle-ductile transition and associated seismicity : Experimental and numerical studies and relationship with the b value, *Journal of Geophysical Research : Solid Earth*, **108**(B1). 17

Ampuero, J.-P., and A. M. Rubin (2008), Earthquake nucleation on rate and state faults-aging and slip laws, *Journal of Geophysical Research : Solid Earth*, **113**(B1). 15, 18

Avouac, J.-P. (2015), From geodetic imaging of seismic and aseismic fault slip to dynamic modeling of the seismic cycle, *Annual Review of Earth and Planetary Sciences*, **43**, 233–271. 10, 12, 15, 16

Bachmann, C. E., S. Wiemer, B. Goertz-Allmann, and J. Woessner (2012), Influence of pore-pressure on the event-size distribution of induced earthquakes, *Geophysical Research Letters*, **39**(9). 17

Bagiya, M. S., D. Thomas, E. Astafyeva, Q. Bletery, P. Lognonné, and D. S. Ramesh (2020), The ionospheric view of the 2011 tohoku-oki earthquake seismic source : the first 60 seconds of the rupture, *Scientific reports*, **10**(1), 1–15. 5

Barbot, S., N. Lapusta, and J.-P. Avouac (2012), Under the hood of the earthquake machine : Toward predictive modeling of the seismic cycle, *Science*, **336**(6082), 707–710. 15

Bartlow, N. M., S. Miyazaki, A. M. Bradley, and P. Segall (2011), Space-time correlation of slip and tremor during the 2009 cascadia slow slip event, *Geophysical Research Letters*, **38**(18). 12

Bergen, K. J., and G. C. Beroza (2018), Detecting earthquakes over a seismic network using single-station similarity measures, *Geophysical Journal International*, **213**(3), 1984–1998. 26

Bergen, K. J., P. A. Johnson, V. Maarten, and G. C. Beroza (2019), Machine learning for data-driven discovery in solid earth geoscience, *Science*, **363**(6433). 24, 26

Bletery, Q. (2015), Analyse probabiliste et multi-données de la source de grands séismes, Ph.D. thesis, Nice. 5, 9
Bletery, Q., and J.-M. Nocquet (2020), Slip bursts during coalescence of slow slip events in cascadia, *Nature communications*, **11**(1), 1–6. 60, 73

Bletery, Q., A. Sladen, B. Delouis, M. Vallée, J.-M. Nocquet, L. Rolland, and J. Jiang (2014), A detailed source model for the mw9.0 tohoku-oki earthquake reconciling geodesy, seismology, and tsunami records, *Journal of Geophysical Research : Solid Earth*, **119**(10), 7636–7653. 8, 9

Bletery, Q., A. Sladen, B. Delouis, and L. Mattéo (2015), Quantification of tsunami bathymetry effect on finite fault slip inversion, *Pure and Applied Geophysics*, **172**(12), 3655–3670. 7

Bletery, Q., A. Sladen, J. Jiang, and M. Simons (2016a), A bayesian source model for the 2004 great sumatra-andaman earthquake, *Journal of Geophysical Research : Solid Earth*, **121**(7), 5116–5135. 9

Bletery, Q., A. M. Thomas, A. W. Rempel, L. Karlstrom, A. Sladen, and L. De Barros (2016b), Mega-earthquakes rupture flat megathrusts, *Science*, **354**(6315), 1027–1031. 28

Bletery, Q., A. M. Thomas, A. W. Rempel, and J. L. Hardebeck (2017a), Imaging shear strength along subduction faults, *Geophysical Research Letters*, **44**(22), 11–329. 35

Bletery, Q., A. M. Thomas, J. C. Hawthorne, R. M. Skarbek, A. W. Rempel, and R. D. Krogstad (2017b), Characteristics of secondary slip fronts associated with slow earthquakes in cascadia, *Earth and Planetary Science Letters*, **463**, 212–220. 62

Bletery, Q., O. Cavalié, J.-M. Nocquet, and T. Ragon (2020), Distribution of interseismic coupling along the north and east anatolian faults inferred from insar and gps data, *Geophysical Research Letters*, **47**(16), e2020GL087,775. 47

Böse, M., E. Hauksson, K. Solanki, H. Kanamori, and T. Heaton (2009), Real-time testing of the on-site warning algorithm in southern california and its performance during the july 29 2008 mw5.4 chino hills earthquake, *Geophysical Research Letters*, **36**(5). 20

Böse, M., T. H. Heaton, and E. Hauksson (2012), Real-time finite fault rupture detector (finder) for large earthquakes, *Geophysical Journal International*, **191**(2), 803–812. 20

Böse, M., C. Felizardo, and T. H. Heaton (2015), Finite-fault rupture detector (finder) : Going real-time in californian shakealert warning system, *Seismological Research Letters*, **86**(6), 1692–1704. 20
Böse, M., D. E. Smith, C. Felizardo, M.-A. Meier, T. H. Heaton, and J. F. Clinton (2018), Finder v. 2 : Improved real-time ground-motion predictions for m2–m9 with seismic finite-source characterization, *Geophysical Journal International*, 212(1), 725–742. 20

Bouchon, M. (1981), A simple method to calculate green's functions for elastic layered media, *Bulletin of the Seismological Society of America*, 71(4), 959–971. 6

Bouchon, M., H. Karabulut, M. Aktar, S. Özalaybey, J. Schmittbuhl, and M.-P. Bouin (2011), Extended nucleation of the 1999 mw 7.6 izmit earthquake, *science*, 331(6019), 877–880. 18

Bouchon, M., V. Durand, D. Marsan, H. Karabulut, and J. Schmittbuhl (2013), The long precursory phase of most large interplate earthquakes, *Nature geoscience*, 6(4), 299–302. 18

Brodsky, E. E., and T. Lay (2014), Recognizing foreshocks from the 1 april 2014 chile earthquake, *Science*, 344(6185), 700–702. 18

Calais, E., and J. B. Minster (1995), Gps detection of ionospheric perturbations following the january 17, 1994, northridge earthquake, *Geophysical Research Letters*, 22(9), 1045–1048. 5

Calais, E., J. Han, C. DeMets, and J. Nocquet (2006), Deformation of the north american plate interior from a decade of continuous gps measurements, *Journal of geophysical research : solid earth*, 111 (B6). 4

Calais, E., S. Symithe, T. Monfret, B. Delouis, A. Lomax, F. Courboulex, J.-P. Ampuero, P. Lara, Q. Bletery, J. Chèze, F. Peix, A. Deschamps, B. de Lépinay, B. Raimbault, R. Jolivet, S. Paul, S. St-Fleur, D. Boisson, Y. Fukushima, Z. Duputel, L. Xu, and L. Meng (2022), Citizen seismology helps decipher the 2021 haiti earthquake, *Science*, 376(6590), 283–287. 114

Cao, T., and K. Aki (1986), Seismicity simulation with a rate-and state-dependent friction law, in *Friction and Faulting*, pp. 487–513, Springer. 27

Cavalié, O., C. Lasserre, M.-P. Doin, G. Peltzer, J. Sun, X. Xu, and Z.-K. Shen (2008), Measurement of interseismic strain across the haiyuan fault (gansu, china), by insar, *Earth and Planetary Science Letters*, 275(3-4), 246–257. 5

Choy, G. L., and S. H. Kirby (2004), Apparent stress, fault maturity and seismic hazard for normal-fault earthquakes at subduction zones, *Geophysical Journal International*, 159(3), 991–1012. 27
Chung, A. I., I. Henson, and R. M. Allen (2019), Optimizing earthquake early warning performance : Elarms-3, *Seismological Research Letters*, 90(2A), 727-743. 20

Crowell, B. W., D. A. Schmidt, P. Bodin, J. E. Vidale, J. Gomberg, J. Renate Hartog, V. C. Kress, T. I. Melbourne, M. Santillan, S. E. Minson, et al. (2016), Demonstration of the cascadia g-fast geodetic earthquake early warning system for the nisqually, washington, earthquake, *Seismological Research Letters*, 87(4), 930-943. 21

Crowell, B. W., D. A. Schmidt, P. Bodin, J. E. Vidale, B. Baker, S. Barrientos, and J. Geng (2018), G-fast earthquake early warning potential for great earthquakes in chile, *Seismological Research Letters*, 89(2A), 542-556. 21

Dai, H., and C. MacBeth (1995), Automatic picking of seismic arrivals in local earthquake data using an artificial neural network, *Geophysical journal international*, 120(3), 758-774. 24

DeMets, C., P. E. Jansma, G. S. Mattioli, T. H. Dixon, F. Farina, R. Bilham, E. Calais, and P. Mann (2000), Gps geodetic constraints on caribbean-north america plate motion, *Geophysical Research Letters*, 27(3), 437-440. 4

Derras, B., P. Y. Bard, and F. Cotton (2014), Towards fully data driven ground-motion prediction models for europe, *Bulletin of Earthquake Engineering*, 12(1), 495-516. 26

DeVries, P. M., T. B. Thompson, and B. J. Meade (2017), Enabling large-scale viscoelastic calculations via neural network acceleration, *Geophysical Research Letters*, 44(6), 2662-2669. 26

Dieterich, J. (1994), A constitutive law for rate of earthquake production and its application to earthquake clustering, *Journal of Geophysical Research : Solid Earth*, 99(B2), 2601-2618. 15

Dieterich, J. H. (1979), Modeling of rock friction : 1. experimental results and constitutive equations, *Journal of Geophysical Research : Solid Earth*, 84(B5), 2161-2168. 14

Dieterich, J. H. (1987), Nucleation and triggering of earthquake slip : effect of periodic stresses, *Tectonophysics*, 144(1-3), 127-139. 15, 18

Doin, M.-P., C. Lasserre, G. Peltzer, O. Cavalié, and C. Doubre (2009), Corrections of stratified tropospheric delays in sar interferometry : Validation with global atmospheric models, *Journal of Applied Geophysics*, 69(1), 35-50. 5
Dowla, F. U., S. R. Taylor, and R. W. Anderson (1990), Seismic discrimination with artificial neural networks : preliminary results with regional spectral data, *Bulletin of the Seismological Society of America*, **80**(5), 1346-1373. 24

Dragert, H., K. Wang, and T. S. James (2001), A silent slip event on the deeper cascadia subduction interface, *Science*, **292**(5521), 1525-1528. 12

Duputel, Z., L. Rivera, H. Kanamori, G. P. Hayes, B. Hirshorn, and S. Weinstein (2011), Real-time w phase inversion during the 2011 off the pacific coast of tohoku earthquake, *Earth, planets and space*, **63**(7), 535-539. 21

Duputel, Z., L. Rivera, H. Kanamori, and G. Hayes (2012), W phase source inversion for moderate to large earthquakes (1990-2010), *Geophysical Journal International*, **189**(2), 1125-1147. 21

Duputel, Z., P. S. Agram, M. Simons, S. E. Minson, and J. L. Beck (2014), Accounting for prediction uncertainty when inferring subsurface fault slip, *Geophysical Journal International*, **197**(1), 464-482. 9

Dysart, P. S., and J. J. Pulli (1990), Regional seismic event classification at the noress array : seismological measurements and the use of trained neural networks, *Bulletin of the Seismological Society of America*, **80**(6B), 1910-1933. 24

Ellsworth, W. L., and G. C. Beroza (1998), Observation of the seismic nucleation phase in the ridgecrest, california, earthquake sequence, *Geophysical research letters*, **25**(3), 401-404. 84

Ellsworth, W. L., and F. Bulut (2018), Nucleation of the 1999 izmit earthquake by a triggered cascade of foreshocks, *Nature Geoscience*, **11**(7), 531-535. 18

Espinosa-Aranda, J. M., A. Cuellar, A. Garcia, G. Ibarrola, R. Islas, S. Maldonado, and F. Rodriguez (2009), Evolution of the mexican seismic alert system (sasmex), *Seismological Research Letters*, **80**(5), 694-706. 83

Frank, W. B., B. Rousset, C. Lasserre, and M. Campillo (2018), Revealing the cluster of slow transients behind a large slow slip event, *Science advances*, **4**(5), eaat0661. 14

Fujinawa, Y., and Y. Noda (2013), Japan's earthquake early warning system on 11 march 2011 : performance, shortcomings, and changes, *Earthquake Spectra*, **29**(1_suppl), 341-368. 83
Geller, R. J. (1997), Earthquake prediction : a critical review, *Geophysical Journal International*, **131**(3), 425-450. 18, 19

Gerstenberger, M. C., S. Wiemer, L. M. Jones, and P. A. Reasenberg (2005), Real-time forecasts of tomorrow's earthquakes in california, *Nature*, **435**(7040), 328-331. 17

Goebel, T., D. Schorlemmer, T. Becker, G. Dresen, and C. Sammis (2013), Acoustic emissions document stress changes over many seismic cycles in stick-slip experiments, *Geophysical Research Letters*, **40**(10), 2049-2054. 17

Gomberg, J., C. 2007, and B. W. Group (2010), Slow-slip phenomena in cascadia from 2007 and beyond : A review, *Bulletin*, **122**(7-8), 963-978. 18

Gomberg, J., A. Wech, K. Creager, K. Obara, and D. Agnew (2016), Reconsidering earthquake scaling, *Geophysical Research Letters*, **43**(12), 6243-6251. 14

Grapenthin, R., I. Johanson, and R. M. Allen (2014a), The 2014 mw 6.0 napa earthquake, california : Observations from real-time gps-enhanced earthquake early warning, *Geophysical Research Letters*, **41**(23), 8269-8276. 21

Grapenthin, R., I. A. Johanson, and R. M. Allen (2014b), Operational real-time gps-enhanced earthquake early warning, *Journal of Geophysical Research : Solid Earth*, **119**(10), 7944-7965. 21

Gulia, L., and S. Wiemer (2010), The influence of tectonic regimes on the earthquake size distribution : A case study for italy, *Geophysical Research Letters*, **37**(10). 17

Gulia, L., and S. Wiemer (2019), Real-time discrimination of earthquake foreshocks and after-shocks, *Nature*, **574**(7777), 193-199. 17, 18

Gulia, L., A. P. Rinaldi, T. Tormann, G. Vannucci, B. Enescu, and S. Wiemer (2018), The effect of a mainshock on the size distribution of the aftershocks, *Geophysical Research Letters*, **45**(24), 13-277. 17

Gutenberg, B., and C. F. Richter (1944), Frequency of earthquakes in california, *Bulletin of the Seismological society of America*, **34**(4), 185-188. 17

Hanks, T. C., and H. Kanamori (1979), A moment magnitude scale, *Journal of Geophysical Research : Solid Earth*, **84**(B5), 2348-2350. 4
Hardebeck, J. L. (2015), Stress orientations in subduction zones and the strength of subduction
megathrust faults, *Science*, **349**(6253), 1213–1216. 59

Harms, J. (2016), Transient gravity perturbations from a double-couple in a homogeneous half-space, *Geophysical Journal International*, **205**(2), 1153–1164. 22

Harms, J., J.-P. Ampuero, M. Barsuglia, E. Chassande-Mottin, J.-P. Montagner, S. Somala, and B. Whiting (2015), Transient gravity perturbations induced by earthquake rupture, *Geophysical Journal International*, **201**(3), 1416–1425. 22

Hayes, G. P., D. J. Wald, and R. L. Johnson (2012), Slab1.0: A three-dimensional model of global subduction zone geometries, *Journal of Geophysical Research: Solid Earth*, **117**(B1). 7

Hayes, G. P., G. L. Moore, D. E. Portner, M. Hearne, H. Flamme, M. Furtney, and G. M. Smoczynyk (2018), Slab2, a comprehensive subduction zone geometry model, *Science*, **362**(6410), 58–61. 7

Heaton, T. H. (1985), A model for a seismic computerized alert network, *Science*, **228**(4702), 987–990. 84

Hecker, S., T. Dawson, and D. Schwartz (2010), Normal-faulting slip maxima and stress-drop variability: A geological perspective, *Bulletin of the Seismological Society of America*, **100**(6), 3130–3147. 28

Heki, K., Y. Otsuka, N. Choosakul, N. Hemmakorn, T. Komolmis, and T. Maruyama (2006), Detection of ruptures of andaman fault segments in the 2004 great sumatra earthquake with coseismic ionospheric disturbances, *Journal of Geophysical Research: Solid Earth*, **111**(B9). 5

Hicks, S. (2015), Seismic properties and processes along the subduction plate interface : the februrary 2010 mw 8.8 maule, chile earthquake, Ph.D. thesis, University of Liverpool. 27

Hino, R., Y. Tanioka, T. Kanazawa, S. Sakai, M. Nishino, and K. Suyehiro (2001), Micro-tsunami from a local interplate earthquake detected by cabled offshore tsunami observation in northeastern japan, *Geophysical Research Letters*, **28**(18), 3533–3536. 5

Holtzman, B. K., A. Paté, J. Paisley, F. Waldhauser, and D. Repetto (2018), Machine learning reveals cyclic changes in seismic source spectra in geysers geothermal field, *Science advances*, **4**(5), eaao2929. 26
Hoshiba, M., and S. Aoki (2015), Numerical shake prediction for earthquake early warning : Data assimilation, real-time shake mapping, and simulation of wave propagation, *Bulletin of the Seismological Society of America*, **105**(3), 1324-1338. 20

Hoshiba, M., K. Iwakiri, N. Hayashimoto, T. Shimoyama, K. Hirano, Y. Yamada, Y. Ishigaki, and H. Kikuta (2011), Outline of the 2011 off the pacific coast of tohoku earthquake (mw 9.0)—earthquake early warning and observed seismic intensity—, *Earth, planets and space*, **63**(7), 547-551. 83

Hulbert, C., B. Rouet-Leduc, P. A. Johnson, C. X. Ren, J. Rivière, D. C. Bolton, and C. Marone (2019), Similarity of fast and slow earthquakes illuminated by machine learning, *Nature Geoscience*, **12**(1), 69-74. 26

Ide, S., G. C. Beroza, D. R. Shelly, and T. Uchide (2007), A scaling law for slow earthquakes, *Nature*, **447**(7140), 76-79. 13, 14, 62

Ito, Y., K. Obara, K. Shiomi, S. Sekine, and H. Hirose (2007), Slow earthquakes coincident with episodic tremors and slow slip events, *Science*, **315**(5811), 503-506. 13

Jiang, J., and N. Lapusta (2016), Deeper penetration of large earthquakes on seismically quiescent faults, *Science*, **352**(6291), 1293-1297. 15

Jolivet, R., and W. Frank (2020), The transient and intermittent nature of slow slip, *AGU Advances*, **1**(1), e2019AV000,126. 13, 14

Jolivet, R., R. Grandin, C. Lasserre, M.-P. Doin, and G. Peltzer (2011), Systematic insar tropospheric phase delay corrections from global meteorological reanalysis data, *Geophysical Research Letters*, **38**(17). 5

Jolivet, R., M. Simons, P. Agram, Z. Duputel, and Z.-K. Shen (2015), Aseismic slip and seismic coupling along the central san andreas fault, *Geophysical Research Letters*, **42**(2), 297-306. 12

Jolivet, R., M. Simons, Z. Duputel, J.-A. Olive, H. Bhat, and Q. Bletery (2020), Interseismic loading of subduction megathrust drives long-term uplift in northern chile, *Geophysical Research Letters*, **47**(8), e2019GL085,377. 12, 60

Jordan, M. I., and T. M. Mitchell (2015), Machine learning : Trends, perspectives, and prospects, *Science*, **349**(6245), 255-260. 24
Jordan, T. H., Y.-T. Chen, P. Gasparini, R. Madariaga, I. Main, W. Marzocchi, G. Papado-
poulos, G. Sobolev, K. Yamaoka, and J. Zschau (2011), Operational earthquake forecasting.
state of knowledge and guidelines for utilization, *Annals of Geophysics*, **54** (4). 19

Juhel, K., J.-P. Ampuero, M. Barsuglia, P. Bernard, E. Chassande-Mottin, D. Fiorucci, J. Harms, J.-P. Montagner, M. Vallée, and B. Whiting (2018), Earthquake early warning using future generation gravity strainmeters, *Journal of Geophysical Research : Solid Earth*, **123**(12), 10–889. 24

Kagan, Y. Y. (1997), Are earthquakes predictable?, *Geophysical Journal International*, **131**(3), 505-525. 18, 19

Kanamori, H. (2005), Real-time seismology and earthquake damage mitigation, *Annu. Rev. Earth Planet. Sci.*, **33**, 195–214. 84

Kanamori, H., and J. W. Given (1981), Use of long-period surface waves for rapid determination of earthquake-source parameters, *Physics of the Earth and Planetary interiors*, **27**(1), 8–31.
21

Kanamori, H., and L. Rivera (2008), Source inversion of wphase : speeding up seismic tsunami warning, *Geophysical Journal International*, **175**(1), 222–238. 21, 22

Kaneko, Y., J.-P. Avouac, and N. Lapusta (2010), Towards inferring earthquake patterns from geodetic observations of interseismic coupling, *Nature Geoscience*, **3**(5), 363–369. 15

Kato, A., K. Obara, T. Igarashi, H. Tsuruoka, S. Nakagawa, and N. Hirata (2012), Propagation of slow slip leading up to the 2011 mw 9.0 tohoku-oki earthquake, *Science*, **335**(6069), 705–708. 18

Käufl, P., A. P. Valentine, and J. Trampert (2016), Probabilistic point source inversion of strong-motion data in 3-d media using pattern recognition : A case study for the 2008 mw 5.4 chino hills earthquake, *Geophysical Research Letters*, **43**(16), 8492–8498. 26

Kawamoto, S., Y. Hiyama, Y. Ohta, and T. Nishimura (2016), First result from the geonet real-time analysis system (regard) : the case of the 2016 kumamoto earthquakes, *Earth, Planets and Space*, **68**(1), 1–12. 21

Kawamoto, S., Y. Ohta, Y. Hiyama, M. Todoriki, T. Nishimura, T. Furuya, Y. Sato, T. Yahagi, and K. Miyagawa (2017), Regard : A new gnss-based real-time finite fault modeling system for geonet, *Journal of Geophysical Research : Solid Earth*, **122**(2), 1324–1349. 21
Kennett, B., and E. Engdahl (1991), *Traveltimes for global earthquake location and phase identification*, *Geophysical Journal International*, **105**(2), 429-465. **7**

Kilb, D., and J. Gomberg (1999), *The initial subevent of the 1994 northridge, california, earthquake : Is earthquake size predictable?*, *Journal of Seismology*, **3**(4), 409-420. **84**

King, R. (1995), *Documentation for the gamit gps analysis software*, *Mass. Inst. of Technol*. **4**

Kohler, M. D., E. S. Cochran, D. Given, S. Guiwits, D. Neuhauser, I. Henson, R. Hartog, P. Bodin, V. Kress, S. Thompson, et al. (2018), *Earthquake early warning shakealert system : West coast wide production prototype*, *Seismological Research Letters*, **89**(1), 99-107. **83**

Kong, Q., R. M. Allen, and L. Schreier (2016a), *Myshake : Initial observations from a global smartphone seismic network*, *Geophysical Research Letters*, **43**(18), 9588-9594. **20**

Kong, Q., R. M. Allen, L. Schreier, and Y.-W. Kwon (2016b), *Myshake : A smartphone seismic network for earthquake early warning and beyond*, *Science advances*, **2**(2), e1501055. **20**, **26**

Lasserre, C., G. Peltzer, F. Crampé, Y. Klinger, J. Van der Woerd, and P. Tapponnier (2005), *Coseismic deformation of the 2001 mw= 7.8 kokoxili earthquake in tibet, measured by synthetic aperture radar interferometry*, *Journal of Geophysical Research : Solid Earth*, **110**(B12). **5**

Lay, T. (2018), *A review of the rupture characteristics of the 2011 tohoku-oki mw 9.1 earthquake*, *Tectonophysics*, **733**, 4-36. **9**

Licciardi, A., Q. Bletery, B. Rouet-Leduc, J.-P. Ampuero, and K. Juhel (2022), *Instantaneous tracking of earthquake growth with elastogravity signals*, *Nature*, **606**, 319-324. **85**

Lin, J.-T., D. Melgar, A. Thomas, and J. Searcy (2021), *Early warning for great earthquakes from characterization of crustal deformation patterns with deep learning*, *Journal of Geophysical Research : Solid Earth*, **126**(10), e2021JB022703. **26**

Loveless, J. P., and B. J. Meade (2011), *Spatial correlation of interseismic coupling and coseismic rupture extent of the 2011 mw= 9.0 tohoku-oki earthquake*, *Geophysical Research Letters*, **38**(17). **12**

Maeda, T., T. Furumura, S. Sakai, and M. Shinohara (2011), *Significant tsunami observed at ocean-bottom pressure gauges during the 2011 off the pacific coast of tohoku earthquake*, *Earth, Planets and Space*, **63**(7), 803-808. **5**
Manighetti, I., M. Campillo, S. Bouley, and F. Cotton (2007), Earthquake scaling, fault segmentation, and structural maturity, *Earth and Planetary Science Letters*, **253**(3-4), 429-438.
**27**, **28**

Manighetti, I., D. Zigone, M. Campillo, and F. Cotton (2009), Self-similarity of the largest-scale segmentation of the faults : Implications for earthquake behavior, *Earth and Planetary Science Letters*, **288**(3-4), 370-381.
**28**

Mansinha, L., and D. Smylie (1971), The displacement fields of inclined faults, *Bulletin of the Seismological Society of America*, **61**(5), 1433-1440.
**6**

Marone, C. (1998), Laboratory-derived friction laws and their application to seismic faulting, *Annual Review of Earth and Planetary Sciences*, **26**(1), 643-696.
**14**

Marone, C. J., C. Scholtz, and R. Bilham (1991), On the mechanics of earthquake afterslip, *Journal of Geophysical Research : Solid Earth*, **96**(B5), 8441-8452.
**15**

Marzocchi, W., M. Taroni, and G. Falcone (2017), Earthquake forecasting during the complex amatrice-norcia seismic sequence, *Science advances*, **3**(9), e1701,239.
**17**

Massonnet, D., M. Rossi, C. Carmona, F. Adragna, G. Peltzer, K. Feigl, and T. Rabaute (1993), The displacement field of the landers earthquake mapped by radar interferometry, *nature*, **364**(6433), 138-142.
**5**

McCann, M. T., K. H. Jin, and M. Unser (2017), Convolutional neural networks for inverse problems in imaging : A review, *IEEE Signal Processing Magazine*, **34**(6), 85-95.
**26**

Meier, M.-A. (2017), How “good” are real-time ground motion predictions from earthquake early warning systems?, *Journal of Geophysical Research : Solid Earth*, **122**(7), 5561-5577.
**21**, **84**

Meier, M.-A., T. Heaton, and J. Clinton (2015), The gutenberg algorithm : Evolutionary bayesian magnitude estimates for earthquake early warning with a filter bank, *Bulletin of the Seismological Society of America*, **105**(5), 2774-2786.
**20**

Meier, M.-A., J. Ampuero, and T. H. Heaton (2017), The hidden simplicity of subduction megathrust earthquakes, *Science*, **357**(6357), 1277-1281.
**23**, **84**

Melgar, D., T. I. Melbourne, B. W. Crowell, J. Geng, W. Szeliga, C. Scrivner, M. Santillan, and D. E. Goldberg (2020), Real-time high-rate gnss displacements : Performance demonstration
during the 2019 ridgecrest, california, earthquakes, *Seismological Research Letters*, 91(4), 1943-1951. 4

Metois, M., C. Vigny, and A. Socquet (2016), Interseismic coupling, megathrust earthquakes and seismic swarms along the chilean subduction zone (38-18 s), *Pure and Applied Geophysics*, 173(5), 1431-1449. 12

Michel, S., A. Gualandi, and J.-P. Avouac (2019), Similar scaling laws for earthquakes and cascadia slow-slip events, *Nature*, 574 (7779), 522-526. 14

Mignan, A. (2015), Modeling aftershocks as a stretched exponential relaxation, *Geophysical Research Letters*, 42(22), 9726-9732. 16

Mignan, A. (2016), Reply to “comment on ‘revisiting the 1894 omori aftershock dataset with the stretched exponential function’ by a. mignan” by s. hainzl and a. christophersen, *Seismological Research Letters*, 87(5), 1134-1137. 16

Miller, M. M., T. Melbourne, D. J. Johnson, and W. Q. Sumner (2002), Periodic slow earthquakes from the cascadia subduction zone, *Science*, 295(5564), 2423-2423. 12

Minson, S., M. Simons, and J. Beck (2013), Bayesian inversion for finite fault earthquake source models i—theory and algorithm, *Geophysical Journal International*, 194(3), 1701-1726. 9

Minson, S. E., J. R. Murray, J. O. Langbein, and J. S. Gomberg (2014), Real-time inversions for finite fault slip models and rupture geometry based on high-rate gps data, *Journal of Geophysical Research : Solid Earth*, 119(4), 3201-3231. 21

Minson, S. E., M.-A. Meier, A. S. Baltay, T. C. Hanks, and E. S. Cochran (2018), The limits of earthquake early warning : Timeliness of ground motion estimates, *Science advances*, 4(3), eaaq0504. 84

Montagner, J.-P., K. Juhel, M. Barsuglia, J. P. Ampuero, E. Chassande-Mottin, J. Harms, B. Whiting, P. Bernard, E. Clévédé, and P. Lognonné (2016), Prompt gravity signal induced by the 2011 tohoku-oki earthquake, *Nature communications*, 7(1), 1-7. 22

Mori, J., and H. Kanamori (1996), Initial rupture of earthquakes in the 1995 ridgecrest, california sequence, *Geophysical Research Letters*, 23(18), 2437-2440. 84

Mousavi, S. M., Y. Sheng, W. Zhu, and G. C. Beroza (2019), Stanford earthquake dataset (stead) : A global data set of seismic signals for ai, *IEEE Access*, 7, 179,464-179,476. 112