# Fast and full characterization of large earthquakes from prompt elastogravity signals

Kévin Juhel, Quentin Bletery, Andrea Licciardi, Martin Vallée, Céline Hourcade, Théodore Michel

► To cite this version:

Kévin Juhel, Quentin Bletery, Andrea Licciardi, Martin Vallée, Céline Hourcade, et al. Fast and full characterization of large earthquakes from prompt elastogravity signals. Communications Earth & Environment, 2024, 5, 10.1038/s43247-024-01725-9. insu-04841455

HAL Id: insu-04841455

https://insu.hal.science/insu-04841455v1

Submitted on 16 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Fast and full characterization of large earthquakes from prompt elastogravity signals

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

Check for updates

Kévin Juhel<sup>1,2</sup>✉, Quentin Bletery<sup>1</sup>, Andrea Licciardi<sup>1</sup>, Martin Vallée<sup>3</sup>, Céline Hourcade<sup>1,2</sup> & Théodore Michel<sup>1,4</sup>

Prompt ElastoGravity Signals are light-speed gravity-induced signals recorded before the arrival of seismic waves. They have raised interest for early warning applications but their weak amplitudes close to the background seismic noise have questioned their actual potential for operational use. A deep-learning model has recently demonstrated its ability to mitigate this noise limitation and to provide in near real-time the earthquake magnitude (*M*<sub>w</sub>). However, this approach was efficient only for large earthquakes (*M*<sub>w</sub> ≥ 8.3) of known focal mechanism. Here we show unprecedented performance in full earthquake characterization using the dense broadband seismic network deployed in Alaska and Western Canada. Our deep-learning model provides accurate magnitude and focal mechanism estimates of *M*<sub>w</sub> ≥ 7.8 earthquakes, 2 minutes after origin time (hence the tsunamigenic potential). Our results represent a major step towards the routine use of prompt elastogravity signals in operational warning systems, and demonstrate its potential for tsunami warning in densely-instrumented areas.

The size of a tsunami is roughly proportional to the volume of elevated water, and therefore primarily depends on the earthquake magnitude and focal mechanism (as some mechanisms are much more efficient at generating uplift than others). Though steep bathymetry may convert horizontal co-seismic sea-floor motion into water uplift<fnref n="1" />, causing strike-slip earthquakes to occasionally generate tsunamis, such a phenomenon is of second-order importance and may be accounted for through pre-computed corrections depending on the source location². Therefore, the critical part of an efficient tsunami warning is the rapid estimation of the magnitude and focal mechanism of the earthquake. Real-time characterization of an event magnitude and focal mechanism is challenging, as seismic waveforms are affected both by the temporal evolution of the moment release and by the spatial complexity of the source. The first seconds of seismic records do not contain enough information to distinguish a moderate magnitude earthquake from a larger one³,<fnref n="4" />, while, for very large events, later wave arrivals typically saturate local broadband sensors and break the point-source approximation. High-rate Global Navigation Satellite Systems (GNSS) measurements are powerful non-saturated observables⁵⁻⁸, but the validity of the point-source approximation has to be carefully evaluated. The use of close GNSS data may require an extended source description, which makes a robust determination more challenging. Thus, for potentially tsunamigenic (*M*<sub>w</sub> ≥ 8) events, reliable estimates of magnitude are difficult to obtain in

near real-time, which results in the inaccurate estimation of the amplitude of the subsequent tsunami wave. Another way to tackle this challenge is to rely on long-period signals recorded at regional distances, correctly modeled with point-source parameters, such as the W-phase⁹ and the recently identified prompt elastogravity signals (coined PEGS). Inversion algorithms relying on the W-phase provide robust solutions down to 10 minutes after the earthquake origin time, using sensors from regional seismological networks¹⁰,¹¹. A PEGS-based source inversion scheme could provide an earlier reliable solution, as PEGS are quickly available after the rupture onset and not prone to saturation.

PEGS are low-frequency signals measurable between the earthquake origin time and the first seismic wave (P-wave) arrival¹²⁻¹⁴. So far, PEGS have been observed on broadband seismometers during six earthquakes with moment magnitude (*M*<sub>w</sub>) ranging from 7.9 to 9.1¹<fnref n="4" />,¹⁵, and modeled with good accuracy by several numerical modeling approaches mostly based on 1D Earth models and point-source approximation¹⁴,¹⁶,¹⁷. Other ground-attached instruments, such as superconducting gravimeters, also have the ability to detect PEGS¹³,¹⁸. Due to their sensitivity to the first-order source parameters (moment magnitude, source location, fault geometry), PEGS can be used for a rapid source determination of large earthquakes¹⁵,¹⁷. However, the very-low signal-to-noise ratios associated with PEGS (whose amplitudes reach a few n.m.s⁻² at most, even for the

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

<footnote n="1">¹Observatoire de la Côte d'Azur, Université Côte d'Azur, IRD, CNRS, Géoazur, Valbonne, France. ²Laboratoire de Planétologie et Géosciences, CNRS UMR 6112, Nantes Université, Université d'Angers, Le Mans Université, Nantes, France. ³Université Paris Cité, Institut de physique du globe de Paris, CNRS, Paris, France.</footnote>

<footnote n="4">⁴Now at Mines Paris, PSL University, CEMEF, CNRS, Paris, France.</footnote>

<footnote>✉e-mail: kjuhel.pro@gmail.com</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 1 | Location of the broadband stations used during training.** Blue squares indicate the location of sensors available during the full deployment of the Alaska USAarray (from 2018 to 2020), while the yellow dots indicate the location of the legacy stations (still available by the end of the 2021 field season). The focal mechanisms show the GCMT solutions<sup>24</sup> for M<sub>w</sub> ≥ 5 earthquakes since 1976 (scaled to moment magnitude, and color-coded based on the focal mechanism type)<sup>54</sup>. The recent large earthquakes studied within this work are explicitly labeled. Black lines show historical fault lines<sup>55</sup>. Cyan lines are the 20 and 30 km isodepths of the subducting slab, according to the Slab2.0 model<sup>56</sup>.</figcaption></figure>

largest events) prevent their routine operational exploitation based on classical source inversion schemes. To our knowledge, a PEGS-based source inversion was only attempted for the 2011 M<sub>w</sub> 9.1 Tohoku-oki earthquake, which had a large final moment magnitude and a fast onset (ensuring an efficient PEGS generation), and benefited from the coverage of good-quality broadband sensors belonging to regional and global seismological networks<sup>17,19,20</sup>.

PEGSNet, a convolutional neural network (CNN) algorithm relying on PEGS recorded by a network of broadband seismometers, has recently been developed for instantaneous magnitude tracking of large earthquakes along the Japanese subduction fault<sup>21</sup>. PEGSNet successfully estimated the moment magnitude accumulated by the 2011 M<sub>w</sub> 9.1 Tohoku-oki earthquake 100 s after its onset. A conservative lower sensitivity to the magnitude was set to M<sub>w</sub> ≥ 8.3 in Japan. PEGSNet was later applied to another tsunamigenic region—the Chilean subduction zone—and estimated a M<sub>w</sub> > 8.7 magnitude for the 2010 M<sub>w</sub> 8.8 Maule earthquake, 90 s after its origin time<sup>22</sup>. Due to limitations associated with the network geometry and the number of available stations, PEGSNet's lower sensitivity was, in that context, limited to M<sub>w</sub> ≥ 8.7.

Besides the earthquake location and magnitude estimations, recent studies reported the application of deep-learning methods to estimate the earthquake source mechanism<sup>23-25</sup>. In this study, we build on PEGSNet and expand its targets to allow the retrieval of the seismic moment tensor. We apply this new algorithm to the Alaska region, a promising case study, as (1) it is a very active area with diverse tectonic settings<sup>26</sup>, and (2) the good existing station coverage has been recently densified by the deployment of the high-quality USAarray network<sup>27,28</sup>.

## CNN-based earthquake determination in Alaska

We generate a dataset of PEGS comprising 512,000 synthetic earthquakes, at the actual locations of the considered broadband sensors. We gather two distinct datasets: a first dataset (hereafter coined “complete”) corresponding to the full deployment of the networks (all good-quality broadband sensors operating between 2018 and 2020), and a second one (hereafter coined “legacy”) corresponding to the legacy sensors (all good-quality broadband

<figure><img src="image_6.png" /><figcaption>**Fig. 2 | Selection of synthetic double-couple mechanisms, extracted from the training set.** The focal mechanisms are scaled to the prescribed moment magnitude, and color-coded based on the focal mechanism type (plunge of tension axis above 45° for thrust faults, plunge of null axis above 45° for strike-slip faults, the remaining sources are labeled as normal faults). The regular grid (dots) shows all possible synthetic source locations.</figcaption></figure>

sensors available from 2018 to 2021, and still operating by the end of 2021). We show in Fig. 1 the locations of the broadband seismometers used in this study, alongside the historical seismicity in the region.

The synthetic sources are pure double-couple mechanisms (strike-slip, thrust, or normal faulting), and their locations are randomly sampled inside geographic regions whose shape is based on the historical seismicity (see Fig. 2 for a selection of sources, and Method and Fig. S1 of the Supplementary Information for more details). We compute the synthetic PEGS waveforms corresponding to this exhaustive set of possible earthquakes using the QSSP algorithm<sup>17,29</sup>. For each synthetic earthquake, we add empirical noise (i.e. real noise recorded by each seismic station) and gather the vertical noisy PEGS waveforms into 2D images (273 or 178 channels, depending on the dataset, × 300 time samples), the last time sample ranging from the synthetic earthquake origin time t<sub>EQ</sub> to t<sub>EQ</sub> + 300 s (see Method and Fig. S2 for details).

PEGSNet's original architecture consists of a sequence of convolutional blocks, followed by fully connected layers<sup>21,22</sup>. We build on PEGSNet's original architecture and expand the dimension of PEGSNet's output layer from three to nine to take into account additional labels (i.e., the six moment tensor components): this allows the model to output a moment tensor solution. Each input image is labeled with the magnitude M<sub>w</sub>(t) of the synthetic earthquake at the time t of the last data sample, and the source latitude and longitude<sup>21,22</sup>, as well as six moment tensor parameters (see Fig. 3 and Method for details).

## Early source estimations for synthetic earthquakes

We use the geometrical similarity α<sup>30</sup> to assess the performance of moment tensor reconstruction. α is a scalar parameter measuring the geometrical difference between two moment tensors, regardless of their scalar moment: the similarity is equal to 0 for opposite moment tensors, and can grow up to 1 in the case of identical moment tensors.

In what follows, we define a magnitude estimation as accurate if the time-dependent moment magnitude M<sub>w</sub><sup>pred</sup>(t) lies within 0.4 magnitude units from the ground truth M<sub>w</sub><sup>true</sup>(t), and a successful mechanism reconstruction if α ≥ 0.8 (compared to the ground truth, input focal mechanism). Note that M<sub>w</sub><sup>true</sup>(t) is not the final magnitude but the time-dependent ongoing magnitude.

Figure 4a–c shows the magnitude regression performances obtained with the complete network, depending on the predicted focal mechanism. We first note that for a given magnitude, the M<sub>w</sub> estimates are more accurate when a normal or a strike-slip earthquake is predicted, compared to a thrust event: the algorithm estimations exceed the 90% accuracy level when the predicted magnitude is above M<sub>w</sub><sup>pred</sup> = 9.0 for thrust, M<sub>w</sub><sup>pred</sup> = 8.2 for
strike-slip and $M_{w}^{pred} = 7.9$ for normal events. This regression behavior can be related to the known property that a subduction thrust earthquake is not an effective source mechanism for PEGS generation¹⁵.

Figure 4d, e shows the corresponding accuracy of mechanism reconstruction. Two minutes after onset time, the algorithm estimations exceed the 75% accuracy level for $M_{w}^{pred} \ge 7.5$ normal earthquakes, $M_{w}^{pred} \ge 7.8$ thrust earthquakes and $M_{w}^{pred} \ge 7.8$ strike-slip earthquakes. For higher magnitudes ($M_{w}^{pred} \ge 8.4$), the reconstruction of predicted strike-slip mechanisms is slightly less accurate than their thrust counterparts. This may be explained by the broader geometric content explored by the strike-slip distributions, compared to the thrust distributions (see Fig. S1). We further illustrate in Figs. S3–S5 the retrieval of the source mechanism (and

location), for synthetic earthquakes whose final magnitude lies within $M_{w} = 8.0 \pm 0.01$, 2 to 4 min after onset time.

We compare in Figs. S6, S7 of the Supplementary Information the regression performances of the complete and legacy network of sensors. We observe similar performances on samples with high signal-to-noise ratios (SNR): the $M_{w}$ regression performance of the complete network is only up to 3% more accurate than the legacy one for very large earthquakes (that is, $M_{w} > 8.9$ for thrust, $M_{w} > 8.2$ for strike-slip and $M_{w} > 7.8$ for normal earthquakes), because the legacy performance is already very good for those events. The higher density of sensors is nonetheless beneficial for the characterization of lower magnitude earthquakes (up to a 10% performance increase in the complete configuration), for which the addition of redundant

<figure><img src="image_7.png" /><figcaption>**Fig. 3** | PEGSNet architecture. The input data is an image of the shape $M \times N$ (purple), where M is the number of time samples and N is the number of sensors. Each convolutional block is composed of a convolutional layer with ReLU activation (orange) and a spatial dropout layer (gray). Max pooling layers (pink) reduce each dimension of the input data by a factor of two. The number of channels used in each</figcaption></figure>

convolutional layer is indicated for clarity. The last convolutional block is connected to dense layers. The output layer uses a tanh activation function to predict 9 values: the moment magnitude ($M_{w}$), the source longitude and latitude ($\phi$ and $\lambda$), and the six components of the moment tensor ($m_{rr}$, $m_{tt}$, ...).

<caption>**Fig. 4** | Accuracy of the source estimations obtained with the complete network of sensors.</caption>

Accuracy of magnitude (a–c) and focal mechanism (d–f) estimations on test set samples. The test set is separated into predicted thrust (a, d), strike-slip (b, e), and normal (c, f) samples. A magnitude estimation is considered successful if $|M_{w}^{true}(t) - M_{w}^{pred}(t)|<0.4$, while a focal mechanism estimation is considered successful if the geometrical similarity $\alpha > 0.8$. Contour lines highlight a 90% (blue), 75% (red), and 60% (green) accuracy. The dashed lines indicate the average $M_{w}^{true}(t)$ for different final magnitudes.

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 5** | Magnitude estimates obtained with the complete network on test set samples. a–d Probability density of the magnitude estimations on samples with true final magnitudes $M_w = 9.0$, $M_w = 7.8$, $M_w = 7.7$ and $M_w = 7.6 \pm 0.05$, respectively. The red lines show the distribution mode (solid) and the 5th and 95th percentiles (dashed). The black lines show the median (solid) and 5th and 95th percentiles (dashed) of the targets, for reference. e Probability density of the magnitudes estimations on all test set events, when PEGS are removed from the synthetic waveforms.</figcaption></figure>

and valuable information from nearby sensors is crucial for an increase of the signal-to-noise ratio. It is even more beneficial for the mechanism reconstruction, as a higher sensor density directly leads to a better resolution on the source radiation pattern. We note that the high level of performance still obtained in the legacy configuration may be explained by an appropriate selection of the removed sensors. Indeed, although the overall sensor density decreases in the legacy configuration - and if we exclude the easternmost part of the network (that is anyway less critical for the detection of thrust and normal events)—both the azimuthal and epicentral coverage of the legacy network remains similar to the complete configuration (see Fig. 1).

To better evaluate the time-dependent performance, we focus on test set samples whose final magnitude lies within the $M_w^{true} = 9.0 \pm 0.05$, $M_w^{true} = 7.8 \pm 0.05$, $M_w^{true} = 7.7 \pm 0.05$ and $M_w^{true} = 7.6 \pm 0.05$ magnitude ranges (Fig. 5a–d for the complete network, see Fig. S8 of the Supplementary Information for the corresponding legacy performances). The algorithm is able to track the ongoing magnitude released by $M_w = 9.0 \pm 0.05$ events in less than 1 min after onset time until the rupture completion. Note that from that time, given the light-speed nature of the messenger signal, the estimation is virtually instantaneous (no time delay between the black and red curves in Fig. 5a). For $M_w = 7.8 \pm 0.05$ events, the model is unable to provide a robust solution during the first 90 s following onset time. From $t = t_{EQ} + 90$ s, the algorithm is able to provide a good estimation of the earthquake's final magnitude: the lower sensitivity of the algorithm to the moment magnitude is thus $M_w = 7.8$ (that is, the minimum magnitude for which the predictions mode reaches the target value).

The distribution mode for $M_w = 7.7 \pm 0.05$ events indeed oscillates between the actual target value and $M_w^{pred} \sim 6.7$. Such a bimodal distribution implies that a set of samples (with low signal-to-noise ratios) is not characterized, while the remaining samples—associated with a low background noise level, an effective source mechanism, and/or an impulsive source-time function (STF)—are accurately determined. We further illustrate the impact of the STF efficiency on our algorithm performance in Fig. S9 of the Supplementary Information. Previous PEGS studies highlighted the crucial role of the STF onset or source duration on PEGS observability<sup>15,17</sup>. A short-duration STF generates larger PEGS than a longer STF, such that an impulsive STF associated with a moderate $M_w$ earthquake may lead to PEGS detection, while a very slow onset—even associated with a very large rupture such as the 2004 $M_w$ 9.3 Sumatra Andaman earthquake—usually compromises a clear detection. We define the STF inefficiency as the ratio of the time needed to release half of the final moment, compared to its corresponding reference functional form<sup>3</sup>. As expected, an STF with a very slow onset tends to hinder the source characterization, while average to fast onsets lead to accurate $M_w$ estimates (for test samples above the sensitivity threshold).

The distribution mode never reaches the target value for $M_w = 7.6 \pm 0.05$ events, that is, well below the algorithm resolving ability (Fig. 5d). The predictions shown in Fig. 5e finally confirm that the estimations rely solely on PEGS-based information: if we remove the PEGS content from the synthetic waveforms (but keep the P-wave arrival information), the model outputs a constant value around $M_w^{pred} = 6.5$. The 99% confidence level of the noise samples lies below the $M_w = 7.8$ sensitivity threshold estimated above at $t \ge t_{EQ} + 60$ s.

## Source parameter retrieval for recent real earthquakes

We test the model on four recent large earthquakes: (1) the July 29, 2021 $M_w$ 8.2 Chignik subduction earthquake<sup>31,32</sup>, (2) the January 23, 2018 $M_w$ 7.9 Kodiak strike-slip earthquake<sup>35,34</sup>, (3) the July 22, 2020 $M_w$ 7.8 Shumagin subduction earthquake<sup>35,36</sup>, and (4) the October 19, 2020 $M_w$ 7.6 Sand Point intraslab earthquake<sup>17,38</sup>. The 2018 $M_w$ 7.9 Kodiak earthquake is the only event for which a PEGS observation has been reported, based on an optimal stacking procedure<sup>15</sup>. The National Tsunami Warning Center (NTWC), in charge of the tsunami monitoring for Alaska, Canada, and the west coast of the United States<sup>39</sup>, issued a warning alert after each of these earthquakes<sup>40,41</sup>. For these events, we process the real waveforms as in the previous section (carefully truncating them at their respective P-wave arrival time, before removing the instrumental response and band-pass filtering), and shape them into 2D images.

In Fig. 6, we compare the obtained solutions to the SCARDEC source-time function database<sup>12</sup> (for the accumulated moment magnitude) and to the GCMT solution (for the focal mechanism). For the 2021 $M_w$ 8.2 Chignik earthquake, the model estimates a $M_w \ge 8.0$ around 60 s after onset time, and a dip-slip mechanism. It converges towards a stabilized set of solutions, starting 120 s after origin time (Fig. 6a), providing faster and more accurate estimations than approaches based on GNSS and seismic data<sup>13,44</sup> (Fig. 6a).

Applied to a significantly smaller thrust rupture (the 2020 $M_w$ 7.8 Shumagin earthquake), the model is still able to provide a good estimate of the event size and a dip-slip mechanism, in less than 140 s (Fig. 6c). This result first makes the Shumagin earthquake the lowest magnitude earthquake for which PEGS has been detected. Second, and more importantly, this detection is associated with an accurate determination of its magnitude and source parameters in very near real-time. This finding demonstrates that the use of PEGS for rapid source characterization is not limited to exceptional events, as $M_w$ 7.8 earthquakes occur on average more than once a year at the global scale. Such a performance is made possible by the large number of stations and by the very efficient learning of the noise-signal separation. It has likely also been helped by the relatively low level of seismic noise observed over the network of sensors prior to the Shumagin earthquake (see Fig. S10).

Despite a relatively high noise level prior to the 2018 $M_w$ 7.9 Kodiak earthquake ($0.36 \text{ nm.s}^{-2}$), the model is also able to estimate an accurate
**Fig. 6 | Source determination for recent large Alaska earthquakes.** Magnitude estimations (in orange) as time grows from onset, for **a** the July 29, 2021 $M_w$ 8.2 Chignik earthquake, **b** the January 23, 2018 $M_w$ 7.9 Kodiak earthquake, **c** the July 22, 2020 $M_w$ 7.8 Shumagin earthquake, and **d** the October 19, 2020 $M_w$ 7.6 Sand Point earthquake. The moment-rate source-time functions from the SCARDEC database¹² are shown in blue, and the corresponding accumulated moment magnitudes are shown in black as a reference for each event. The estimated moment tensor solutions are represented by the orange focal mechanism plots, and the GCMT solution is shown in black as a reference for each event. GNSS-based $M_w$ estimates are also shown in green for comparison (a, b)¹⁴,¹⁵. The top panel relies on the legacy channels, while the bottom three panels rely on the complete network. For each sub-plot, the gray vertical lines indicate P-wave arrival times at the 10th and 100th stations closest to the source.

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

moment magnitude and strike-slip rupture geometry (Fig. 6b), starting 120 s after onset time. As for the 2021 $M_w$ 8.2 Chignik earthquake, the deep-learning solution provides faster and more accurate magnitude estimates than a GNSS-based approach¹⁵.

Finally, the 2020 $M_w$ 7.6 Sand Point earthquake appears to be below the model's lower sensitivity, such that the model is unable to provide an accurate and timely solution. We note that besides the significantly smaller final seismic moment, the 2020 $M_w$ 7.6 Sand Point earthquake also had a slower onset: 15 s after onset time, the event only reached $M_w$ = 7.0 (see the SCARDEC reference curves of Fig. 6), which dampens the efficiency of PEGS generation and makes the accurate retrieval of the source parameters even more challenging.

We further corroborate our findings by generating synthetic PEGS according to the GCMT parameters of the real earthquakes described in this section, and corrupting the waveforms with 1000 time-windows of empirical noise randomly extracted from the test set. Gathering the magnitude predictions for each of these newly generated test samples provides uncertainty estimates for each earthquake as a function of time (see Fig. S11).

We finally assess the benefit of using a deep-learning approach for robust and fast characterization of these large earthquakes, compared to a classical source inversion scheme. We reproduce the PEGS-based, linear inversion documented for the 2011 $M_w$ 9.1 Tohoku earthquake²⁰, and conduct inversions for the 2018 $M_w$ 7.9 Kodiak and 2021 $M_w$ 8.2 Chignik earthquakes. To evaluate uncertainties, we apply both approaches to 1000 test samples of PEGS synthetics augmented with 1000 windows of empirical noise. We find that the source inversions based on PEGS-alone are not accurate enough to constrain the source characteristics (Fig. 7). The hindering effect of seismic noise, apprehended and dampened in the deep-learning approach, indeed leads to an overestimation of the magnitude and unreliable focal mechanisms in the linear inversion scheme, highlighting the need to combine a deep-learning method to a PEGS-based approach to achieve the performance we obtain in Fig. 6.

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

**Fig. 7 | Comparison of PEGS-based source estimations obtained by linear inversion or deep-learning.** Linear inversion (red) and deep-learning (green) magnitude and moment tensor solutions, for **a** the 2021 $M_w$ 8.2 Chignik earthquake and **b** the 2018 $M_w$ 7.9 earthquake. The thick colored lines and focal mechanisms correspond to solutions obtained using the real recorded waveforms, while the density distributions illustrate the solutions behavior for 1000 synthetic samples, obtained by the addition of PEGS corresponding to the earthquakes (based on their GCMT parameters) and 1000 noise time-windows extracted from the test set. The dashed lines show the 25th and 95th percentiles. The accumulated moment magnitude obtained from the SCARDEC database¹² and the GCMT solution are shown in black as a reference for each event.

## Conclusions

We built on the deep-learning algorithm PEGSNet to estimate the magnitude and—for the first time—the focal mechanism of large earthquakes in near real-
time from PEGS, and tested the algorithm in the densely-instrumented region of Alaska. The model outputs reliable solutions down to magnitudes key to mitigate risks associated with tsunamis: the algorithm is able to estimate the magnitude of synthetic earthquakes down to $M_w = 7.8$, 2 minutes after their onset time, and to characterize their focal mechanism (thrust, strike-slip or normal faulting). The model is able to provide robust solutions for three earthquakes that recently occurred in the region, in the $M_w$ [7.8-8.2] magnitude range. We emphasize that there is no upper bound in magnitude that would require a modification of the PEGS-based approach. In contrast, dealing with a wide range of magnitudes with seismic waves usually requires special care or an adaptation of the method (either in the selection of the set of sensors, to ensure a valid point-source approximation, or in the data themselves, to avoid saturation). Our method could thus contribute to the ongoing efforts to improve and refine tsunami warnings in the region$^{46,47}$. The performance strongly relies on the dense network of broadband sensors recently deployed in the region. These performances could be achieved in other earthquake-prone areas, and provide an additional motivation to further instrument these regions.

## Method

As actual PEGS observations are limited in number, we compute a large and realistic synthetic PEGS database in order to train the convolutional neural network. To do so, we (1) compute synthetic PEGS waveforms generated by an exhaustive set of possible earthquake sources, (2) download and process an empirical noise database recorded by the considered set of broadband sensors, and (3) gather the obtained waveforms into 2D labeled images.

## Synthetic earthquake simulation, based on historical seismicity

To simulate the synthetic earthquakes, we first retrieve the historical shallow seismicity (<60 km) reported by Global CMT in Alaska since 1976. We separate the resulting database into strike-slip faulting (plunge of null axis above 45°), thrust faulting (plunge of tension axis above 45°), the remaining sources are labeled as normal faults. We focus on the three geographic zones described by the colored polygons in Fig. 2 of the main text, which contain the majority of the strike-slip (in red), thrust interface (in dark cyan), and extensional outer rise seismicity. The density distributions of the corresponding strike, dip, and rake angles are shown in Fig. S1 of the Supplementary Information. Based on these empirical distributions, we set the characteristics of a 'realistic' synthetic seismicity:

1. To emulate a strike-slip seismicity, we extract the strike angle from a uniform distribution $U(0°, 360°)$, the dip angle from a triangular distribution (lower limit = 55°, upper limit and mode = 90°), and the rake angle from two normal distributions $N(μ = 0°/180°, σ = 15°)$.

2. To emulate a thrust seismicity representative of the Alaska-Aleutian subduction zone, we extract the strike angle from a normal distribution $N(μ = 250°, σ = 15°)$, the dip angle from a normal distribution $N(μ = 25°, σ = 7°)$, and the rake angle from a normal distribution $N(μ = 90°, σ = 20°)$.

3. To emulate a normal seismicity representative of the extensional context in the Alaska-Aleutian subduction fore-arc, we extract the strike angle from a normal distribution $N(μ = 250°, σ = 15°)$, the dip angle from a normal distribution $N(μ = 45°, σ = 7°)$, and the rake angle from a normal distribution $N(μ = 270°, σ = 10°)$.

## Source-time function (STF) database

We compute the moment-rate source-time function (STF) database using the functional form $y(t) = μt \exp(-1/2(\lambda t)^2)$ (where the initial slope $μ$ and the characteristic time scale $λ(M_w)$ are optimized from empirical STF datasets)$^3$. The final moment magnitude $M_w$ is extracted from three uniform distributions: $U(5.5, 8.4)$ for normal, $U(5.5, 8.7)$ for strike-slip and $U(5.5, 9.5)$ for thrust events. We account for the variability of STF duration for a given final magnitude, and add a multiplicative Brownian noise to replicate observed STF fluctuations. A large selection of source-time functions from the training set is shown in Fig. S12.

## Data selection, and processing of noise and PEGS waveforms

The deployment of the Alaska Transportable Array (code TA) was completed in 2017, and the network was operated by the Incorporated Research Institutions for Seismology (IRIS) from 2018 through 2020. It included the installation of 194 new stations and the upgrade of 32 existing broadband stations, leading to a regular grid of broadband seismometers with an 85 km station spacing$^4$. Around 110 TA stations were transitioned to other networks in 2019 and 2020 as part of a legacy network, mostly to the Alaska Earthquake Center (code AK) and the Alaska Volcano Observatory (code AV), while the remaining stations were removed during the 2021 field season.

We retrieve all the broadband data publicly available from 2018 to 2021, inside the geographic region shown in Fig. 1. Data processing is performed as follows$^{14}$: we remove the linear trend and the instrumental response from day-long vertical records (BHZ channels), band-pass filter the resulting traces between 2.0 mHz (high-pass Butterworth causal filter, with two poles) and 30.0 mHz (low-pass Butterworth causal filter, with six poles), and finally decimate the traces down to 1 Hz. Since PEGS are very-low amplitude signals, we remove from the dataset the noisiest sensors, that is, we only keep sensors whose median hourly standard deviation does not exceed $1 \text{ nm.s}^{-2}$.

We use the QSSP algorithm$^{17,29}$ to compute the PEGS Green's functions in between each source and sensor locations, inside a 1D Earth model (AK135). After the convolution with the prescribed source-time function, the PEGS synthetic waveforms are band-pass filtered between 2.0 and 30.0 mHz.

## Architecture, labeling, and training

We build on PEGSNet's original architecture$^{21,22}$, which consists of a sequence of convolutional blocks, followed by fully connected layers (detailed schematic of the architecture in Fig. 3). Each convolutional block is composed of a 2D convolutional layer (with a kernel of size = 3, stride = 1, and padding = 1), a ReLU activation layer and a 2D dropout layer (with probability $p = 0.04$). The successive convolutional layers are equivalent to a series of filters, whose trainable weights are optimized during the CNN training to extract the most meaningful features.

To get the final data waveform, we extract a random noise time window from the empirical noise database, and add it to the PEGS waveforms. We set each trace to zero after their respective P-wave arrival, clip the waveforms to a threshold value ($\pm 10 \text{ nm.s}^{-2}$) in order to limit the influence of very noisy traces, and finally normalize the resulting traces by the threshold value to ensure an easier convergence of the optimizer. We gather the vertical noisy PEGS waveforms into 2D images (273 or 178 channels × 300 time samples), the last time sample ranging from the synthetic earthquake origin time $t_{EQ}$ to $t_{EQ} + 300$ s. We order the waveforms such that close-by sensors have relatively close ranks inside the 2D image (see Fig. S2 of the Supplementary Information).

Each input image is labeled with the time-dependent magnitude $M_w(t)$ (corresponding to the time $t$ of the last sample in the input data sample) and the source latitude and longitude$^{21,22}$. We expand the dimension of the output layer of PEGSNet from three to nine to take into account additional labels (i.e., the six moment tensor components): this allows the model to output a moment tensor solution. Each label is later normalized by its respective extreme values during training.

We split the synthetic PEGS database into 70% training, 20% validation, and 10% test sets. The weights of the neural networks are optimized during training (throughout 200 epochs, using the Adam optimizer) by minimizing the smooth L1 loss (the cutoff for $β = 0.3$) between the targets and the output values. The loss curves corresponding to the training of the neural networks shown within this study are shown in Fig. S13.

## Data availability

The broadband seismic data from the following networks were used in this study (alphabetic order): the Alaska Geophysical Network (AK, https://doi.org/10.7914/SN/AK), the National Tsunami Warning Center Alaska Seismic Network (AT, 10.7914/SN/AT), the Alaska Volcano Observatory (AV,
10.7914/SN/AV), the Canadian National Seismograph Network (CN, 10.7914/SN/CN), the Global Seismograph Network (II, 10.7914/SN/II), the Global Seismograph Network (IU, 10.7914/SN/IU), the USArray Transportable Array (TA, 10.7914/SN/TA). Waveforms were retrieved from IRIS using the Python toolbox ObspyDMT¹⁸, and processed with the Python toolbox ObsPy⁴⁹. SCARDEC source-time functions were accessed on http://scardec projects.sismo.ipgp.fr.

## Code availability

PEGSNet⁵⁰ is built and trained using PyTorch⁵¹. The Python scripts used to generate the input PEGS database and train the neural networks described within this paper are available at the following GitLab repository (https://gitlab.com/kjuhel/PEGSNet). Cartopy (https://scitools.org.uk/cartopy), Matplotlib⁵², and perceptually uniform scientific color maps⁵³ were used for plotting purposes.

Received: 25 March 2024; Accepted: 24 September 2024;

Published online: 04 October 2024

## References

1. Tanioka, Y. & Satake, K. Tsunami generation by horizontal displacement of ocean bottom. *Geophys. Res. Lett.* **23**, 861–864 (1996).

2. Bletery, Q., Sladen, A., Delouis, B. & Mattéo, L. Quantification of tsunami bathymetry effect on finite fault slip inversion. *Pure Appl. Geophys.* **172**, 3655–3670 (2015).

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

4. Renou, J., Vallée, M. & Dublanchet, P. How does seismic rupture accelerate? observational insights from earthquake source time functions. *J. Geophys. Res. Solid Earth* **124**, 8942–8952 (2019).

5. 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).

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

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

8. 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).

9. Kanamori, H. W phase. *Geophys. Res. Lett.* **20**, 1691–1694 (1993).

10. Kanamori, H. & Rivera, L. Source inversion of W phase: speeding up seismic tsunami warning. *Geophys. J. Int.* **175**, 222–238 (2008).

11. Duputel, Z. et al. Real-time W phase inversion during the 2011 Tohoku-oki earthquake. *Earth Planets Space* **63**, 535–539 (2011).

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

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

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

15. 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).

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

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

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

19. Kimura, M. et al. Determination of the source parameters of the 2011 Tohoku-Oki earthquake from three-component pre-P gravity signals recorded by dense arrays in Japan. *Earth Planets Space* **73**, 223 (2021).

20. Juhel, K., Duputel, Z., Rivera, L. & Vallée, M. Early source characterization of large earthquakes using W phase and prompt elastogravity signals. *Seismol. Res. Lett.* **95**, 1558–1568 (Seismological Society of America, 2024).

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

22. Arias, G. et al. Rapid source characterization of the Maule earthquake using prompt elastogravity signals. *J. Geophys. Res. Solid Earth* **128**, e2023JB026588 (2023).

23. Kuang, W., Yuan, C. & Zhang, J. Real-time determination of earthquake focal mechanism via deep learning. *Nat. Commun.* **12**, 1432 (2021).

24. Steinberg, A., Vasyura-Bathke, H., Gaebler, P., Ohrnberger, M. & Ceranna, L. Estimation of seismic moment tensors using variational inference machine learning. *J. Geophys. Res. Solid Earth* **126**, e2021JB022685 (2021).

25. Zhang, H., Innanen, K. A. & Eaton, D. W. Inversion for shear-tensile focal mechanisms using an unsupervised physics-guided neural network. *Seismol. Soc. Am.* **92**, 2282–2294 (2021).

26. Doser, D. I. & Lomas, R. The transition from strike-slip to oblique subduction in southeastern Alaska from seismological studies. *Tectonophysics* **316**, 45–65 (2000).

27. Busby, R. W. & Aderhold, K. The Alaska transportable array: as built. *Seismol. Res. Lett.* **91**, 3017–3027 (2020).

28. Ruppert, N. A. & West, M. E. The impact of USArray on earthquake monitoring in Alaska. *Seismol. Res. Lett.* **91**, 601–610 (2020).

29. Wang, R., Heimann, S., Zhang, Y., Wang, H. & Dahm, T. Complete synthetic seismograms based on a spherical self-gravitating Earth model with an atmosphere-ocean-mantle-core structure. *Geophys. J. Int.* **210**, 1739–1764 (2017).

30. Rivera, L. & Kanamori, H. Diagnosing source geometrical complexity of large earthquakes. *Pure Appl. Geophys.* **171**, 2819–2840 (2014).

31. Ye, L. et al. Rupture model for the 29 July 2021 MW 8.2 Chignik, Alaska earthquake constrained by seismic, geodetic, and tsunami observations. *J. Geophys. Res. Solid Earth* **127**, e2021JB023676 (2022).

32. Liu, C., Lay, T. & Xiong, X. The 29 July 2021 MW 8.2 Chignik, Alaska peninsula earthquake rupture inferred from seismic and geodetic observations: re-rupture of the western 2/3 of the 1938 rupture zone. *Geophys. Res. Lett.* **49**, e2021GL096004 (2022).

33. Ruppert, N. A. et al. Complex faulting and triggered rupture during the 2018 MW 7.9 offshore Kodiak, Alaska, earthquake. *Geophys. Res. Lett.* **45**, 7533–7541 (2018).

34. Krabbenhoeft, A., von Huene, R., Miller, J. J., Lange, D. & Vera, F. Strike-slip 23 January 2018 MW 7.9 Gulf of Alaska rare intraplate earthquake: complex rupture of a fracture zone system. *Sci. Rep.* **8**, 13706 (2018).

35. Crowell, B. W. & Melgar, D. Slipping the Shumagin gap: a kinematic coseismic and early afterslip model of the Mw 7.8 Simeonof Island, Alaska, earthquake. *Geophys. Res. Lett.* **47**, e2020GL090308 (2020).

36. Liu, C., Lay, T., Xiong, X. & Wen, Y. Rupture of the 2020 MW 7.8 earthquake in the Shumagin gap inferred from seismic and geodetic observations. *Geophys. Res. Lett.* **47**, e2020GL090806 (2020).

37. Herman, M. W. & Furlong, K. P. Triggering an unexpected earthquake in an uncoupled subduction zone. *Sci. Adv.* **7**, eabf7590 (2021).

38. Santellanes, S. R., Melgar, D., Crowell, B. W. & Lin, J.-T. Potential megathrust co-seismic slip during the 2020 Sand Point, Alaska strike-slip earthquake. *ESS Open Archive* (2022).
39. Whitmore, P. et al. NOAA/West coast and Alaska tsunami warning center Pacific Ocean response criteria. *Sci. Tsunami Hazards* **27**, 1–19 (2008).

40. Suleimani, E. & West, M. E. Tsunami of the 2018 Mw 7.9 Kodiak Earthquake. In *AGU Fall Meeting Abstracts*, vol. 2019, NH43E-0979 (2019).

41. Becker, N. et al. The Pacific Tsunami Warning Center's response to the 2021 Chignik, Alaska earthquake and its tsunami. In *AGU Fall Meeting Abstracts*, vol. 2021, S55G-0238 (2021).

42. Vallée, M. & Douet, V. A new database of source time functions (STFs) extracted from the SCARDEC method. *Phys. Earth Planet. Inter.* **257**, 149–157 (2016).

43. Parameswaran, R. M., Grapenthin, R., West, M. E. & Fozkos, A. Interchangeable use of GNSS and seismic data for rapid earthquake characterization: 2021 Chignik, Alaska, earthquake. *Seismol. Soc. Am.* **94**, 1367–1378 (2023).

44. Fang, R. et al. Earthquake magnitude scaling using peak ground velocity derived from high-rate GNSS observations. *Seismol. Soc. Am.* **92**, 227–237 (2021).

45. Hodgkinson, K. M., Mencin, D. J., Feaux, K., Sievers, C. & Mattioli, G. S. Evaluation of earthquake magnitude estimation and event detection thresholds for real-time GNSS networks: examples from recent events captured by the network of the Americas. *Seismol. Res. Lett.* **91**, 1628–1645 (2020).

46. Nicolsky, D. J., Suleimani, E. & Koehler, R. D. *Tsunami Inundation Maps for the City of Sand Point, Alaska* (State of Alaska, Department of Natural Resources, Division of Geological ..., 2017).

47. Williamson, A. & Allen, R. M. Improving efficacy of tsunami warnings along the west coast of the United States. *Pure Appl. Geophys.* **180**, 1661–1678 (2023).

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

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

50. Juhel, K., Licciardi, A. & Bletery, Q. PEGSNet: CNN-based rapid earthquake characterization using prompt elastogravity signals (2024).

51. Paszke, A. et al. Pytorch: an imperative style, high-performance deep learning library. In *Proc. 33rd International Conference on Advances in Neural Information Processing Systems* 8026–8037 (Curran Associates Inc., 2019).

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

53. Crameri, F. Scientific colour maps (8.0.1) (2023).

54. Ekström, G., Nettles, M. & Dziewoński, A. The global CMT project 2004–2010: centroid-moment tensors for 13,017 earthquakes. *Phys. Earth Planet. Inter.* **200**, 1–9 (2012).

55. Koehler, R., Burns, R.-E. & PAC Combellick, R. Quaternary faults and folds in Alaska: a digital database (2012).

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

## Acknowledgements

We thank the editor and three anonymous reviewers for their constructive comments and suggestions. We thank Rongjiang Wang for his open-source

QSSP algorithm used to compute PEGS synthetics. We thank Kathleen M. Hodgkinson and Revathy M. Parameswaran for sharing their GNSS-based $M_w$ estimates. This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreements 949221). This work was granted access to the HPC resources of IDRIS under the allocations AP011012126, AD011012142R1, A0121012314, A0101012314, AP011012536, and AD011012142 made by GENCI.

## Author contributions

K.J. and Q.B. had the original idea. K.J. performed the CNN training, made the figures and wrote the core of the article. Q.B. supervised the study, A.L. designed the original CNN and M.V. the synthetic seismicity, C.H. performed alternate trainings, and T.M. optimized the algorithm. All the authors contributed to the writing of 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/s43247-024-01725-9.

**Correspondence** and requests for materials should be addressed to Kévin Juhel.

**Peer review information** *Communications Earth & Environment* thanks Shenjian Zhang, He Tang and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Luca Dal Zilio and Carolina Ortiz Guerrero. A peer review file is available.

**Reprints and permissions 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.

**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, 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 licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence 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 licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2024