# Magnitude estimation and ground motion prediction to harness fiber optic distributed acoustic sensing for earthquake early warning

Itzhak Lior, Diane Rivet, Jean Paul Ampuero, Anthony Sladen, Sergio Barrientos, Rodrigo Sánchez-Olavarría, German Alberto Villarroel Opazo, Jose Antonio Bustamante Prado

► To cite this version:

Itzhak Lior, Diane Rivet, Jean Paul Ampuero, Anthony Sladen, Sergio Barrientos, et al. Magnitude estimation and ground motion prediction to harness fiber optic distributed acoustic sensing for earthquake early warning. Scientific Reports, 2023, 13 (1), pp.424. 10.1038/s41598-023-27444-3. hal-03977430

HAL Id: hal-03977430

https://hal.science/hal-03977430

Submitted on 22 May 2023

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>

Check for updates

OPEN

# Magnitude estimation and ground motion prediction to harness fiber optic distributed acoustic sensing for earthquake early warning

Itzhak Lior¹✉, Diane Rivet², Jean-Paul Ampuero², Anthony Sladen², Sergio Barrientos³, Rodrigo Sánchez-Olavarría³, German Alberto Villarroel Opazo⁴ & Jose Antonio Bustamante Prado⁴

Earthquake early warning (EEW) systems provide seconds to tens of seconds of warning time before potentially-damaging ground motions are felt. For optimal warning times, seismic sensors should be installed as close as possible to expected earthquake sources. However, while the most hazardous earthquakes on Earth occur underwater, most seismological stations are located on-land; precious seconds may go by before these earthquakes are detected. In this work, we harness available optical fiber infrastructure for EEW using the novel approach of distributed acoustic sensing (DAS). DAS strain measurements of earthquakes from different regions are converted to ground motions using a real-time slant-stack approach, magnitudes are estimated using a theoretical earthquake source model, and ground shaking intensities are predicted via ground motion prediction equations. The results demonstrate the potential of DAS-based EEW and the significant time-gains that can be achieved compared to the use of standard sensors, in particular for offshore earthquakes.

While earthquake prediction remains out of reach, continuous seismic monitoring has enabled earthquake early warning (EEW) systems that provide alerts to population centers and critical infrastructure seconds to tens of seconds before intense ground shaking is felt¹⁻⁴. Following rupture initiation, warning may be issued by analyzing recorded ground motions in real-time to assess the earthquake's damage potential. The performance of EEW systems largely depends on the spatial distribution of available seismic sensors⁵; for fast and robust warning issuance, seismic instruments should be installed in proximity to active faults, where earthquakes are expected to occur. While most of the largest and most hazardous earthquakes on Earth occur offshore in subduction zones, the vast majority of seismic stations are located on-land. Thus, valuable time may be lost waiting for seismic waves to reach on-land stations⁶. Current solutions, such as densifying on-land seismic networks and installing cabled ocean bottom sensor networks, are implemented in Japan⁷ and Canada⁸. However, high costs preclude their worldwide implementation. An alternative is to convert existing fiber optic cables into dense seismic networks via the novel technology of distributed acoustic sensing (DAS)⁹,¹⁰. The ever-growing worldwide deployment of optical fiber telecommunication infrastructure, in particular submarine cables, opens opportunities for widespread low-cost implementation of DAS for EEW, circumventing costly ocean-bottom deployments and operations. The potential of seafloor DAS for EEW has not been quantitatively demonstrated yet, a gap addressed in this work.

Over the past several years, the unique advantages of DAS have proven valuable for various seismological purposes including earthquake analysis¹¹⁻¹³ and subsurface imaging¹⁴⁻¹⁷. DAS enables the measurement of transient ground deformations (strains or strain-rates) along optical fibers (e.g., internet cables), such as those currently traversing most of our planet, both on-land and underwater. Unlike point-sensors (e.g., seismometers, accelerometers, GNSS), DAS yields spatially dense longitudinal deformation measurements (every several meters, typically 10) along tens-of-kilometers long optical fibers with a maximum range between 80 and 150 km, depending on the specific DAS interrogator. This technology allows for continuous monitoring of large regions and provides a more complete picture of the seismic wave-field. Measurements are obtained using an interrogator unit which is placed at one end of the cable and sends laser pulses along the fiber. Due to small heterogeneities within the fiber, a fraction of the transmitted light is backscattered via Rayleigh scattering. When seismic waves

¹Institute of Earth Sciences, The Hebrew University, Jerusalem, Israel. ²Université Côte d'Azur, Observatoire de La Côte d'Azur, CNRS, IRD, Géoazur, Valbonne, France. ³Centro Sismológico Nacional, Universidad de Chile, Santiago, Chile. ⁴Gtd Grupo S.A., Santiago, Chile. ✉email: itzhak.lior@mail.huji.ac.il
perturb the cable, heterogeneities along the fiber change positions, and so does the Rayleigh backscattering pattern. The backscattering phase differences between time samples are then translated into strain or strain-rate measurements at spacings of several meters along tens-of-kilometers long fibers¹⁸. This technique allows for the transformation of any optical fiber into a dense array of seismo-acoustic sensors, producing measurements with unprecedented spatial and temporal resolutions.

DAS has key features that are ideally suited for the challenges of EEW. It facilitates spatially and temporally continuous seismic measurements at hard-to-reach places, such as underwater¹⁹ and in boreholes²⁰, closer to earthquake hypocenters. The dense spatial sampling facilitates more reliable separation between earthquakes and noise²¹ compared to point-sensors. Furthermore, the DAS interrogator is sensing the whole fiber from one of its ends, nullifying power and telemetry considerations to distant fiber segments. Thus, the use of optical fibers as dense seismic networks could be decisive in the performance of EEW systems, significantly improving earthquake warning times and allowing for better preparedness for intense shaking.

While the advantages of DAS for EEW are appealing, this novel technology suffers several drawbacks that need to be addressed. For DAS strain measurements to represent ground deformations, fibers should be adequately coupled to the ground. However, since many fibers were not deployed with seismological applications in mind, coupling can vary, and be insufficient along specific sections, for reliable measurements. State-of-the-art DAS interrogator units can sense fibers of up to ~150 km (as demonstrated by the earthquake recording performed by an Alcatel OptoDAS interrogator unit in Supplementary Fig. 1) or up to the first repeater, such that more than one system and fiber may be needed to cover vast regions. DAS records strains or strain-rates: these measurements are very sensitive to the local velocity structure beneath the fiber¹³ and to lateral subsurface heterogeneities¹²,²². This property is troublesome to both earthquake magnitude estimation, which typically relies on ground motion (i.e., displacements, velocities or accelerations) measurements²³, and for earthquake location efforts, where sensitivity to the local subsurface structure may complicate earthquake sources' locations¹². In addition, DAS measures strains along the fiber's axis, such that P-waves recorded by horizontal fibers typically induce low amplitudes that may even be below instrumental noise levels, potentially hindering earthquake location capabilities. Furthermore, while EEW is intended for large earthquakes at short distances, where damages are expected to occur, such DAS observations are currently unavailable. Thus, the effects of DAS amplitude saturation and cable-ground coupling behavior during intense deformations are insufficiently reported and studied. In the following sections we tackle several of the mentioned disadvantages, yet additional work is needed to address all issues to reliably use DAS for EEW.

In this work, we propose the first quantitative real-time schemes that will be part of an operational DAS-based EEW system. Early warning is typically achieved by (1) detecting an earthquake, (2) determining its location, (3) resolving the earthquake source parameters (magnitude and stress drop²⁴), and (4) predicting ground shaking intensities, typically peak ground velocities (PGV) and peak ground accelerations (PGA)²⁵. To the best of our knowledge, these four real-time objectives are yet to be addressed using DAS. Real-time earthquake detection and location may be achieved using either well-established point-sensor-based approaches²⁶⁻²⁹ applied to single or multiple DAS channels, or array processing techniques such as beamforming¹²,³⁰,³¹. While detection can be achieved with relative ease even with point-sensor-based algorithms²⁷,²⁹, earthquake location poses several challenges that are unique to DAS data¹²,²². The recorded strain wavefield may not be coherent enough for reliable earthquake location or may be dominated by scattered waves. In addition, the geographical locations of measurements along the fiber need to be calibrated to reduce earthquake location errors. These issues will be considered when devising real-time earthquake location schemes, a subject of subsequent manuscripts. Here, we address the last two objectives: real-time magnitude estimation and shaking intensity prediction.

Most operational EEW systems rely on empirical relations for both magnitude estimation and ground motion prediction³²,³³. The robustness of these relations largely relies on the quality, quantity, and magnitude range of available earthquake observations³⁴. Since DAS is a relatively new seismic measurement technology⁹, current earthquake DAS datasets are insufficient in all aspects to devise robust empirical methods, and a physics-based approach that does not rely on data availability should be developed²⁴,³⁵,³⁶. Recently, a holistic physics-based approach for real-time earthquake source parameter (magnitude and stress drop) estimation and ground motion prediction has been proposed²⁴. A similar method, adapted to DAS data, is developed here by deriving a theoretical expression for real-time magnitude estimation using the root-mean-squares (rms) of ground accelerations.

Since DAS measures strains (or strain-rates) and earthquake magnitude is directly related to ground motions (displacements, velocities and accelerations)²³, DAS measurements first need to be converted to ground motions³⁷. This objective is typically achieved by using the apparent slowness (reciprocal of velocity), $p_x$, measured along the fiber³⁸:

$$ \frac{d^n}{dt^n} D(t) = \frac{d^{n-1}}{dt^{n-1}} \epsilon(t) / p_x' \quad (1) $$

where $\epsilon(t)$ and $D(t)$ are the time-series of strains and ground displacements, respectively, and $n$ equals 0, 1 or 2 corresponds to conversions to ground displacements, velocities or accelerations, respectively. This relation assumes perfect coupling between the fiber and the Earth, an assumption that was found to hold well for different fibers and installations¹¹,¹³,³⁷,³⁷⁻³⁹. Slowness has been observed to change rapidly both in time and space (along the fiber): temporal variations are due to velocity differences among recorded seismic phases (i.e., P-, S-, surface-waves)³⁷ and spatial variations are a result of lateral subsurface-velocity heterogeneities, that may be significant and abrupt¹²,¹⁶,²²,³⁷. Accurate conversion of strain-rates to ground accelerations requires that slowness be resolved as a function of both time and space. Recently, a slant-stack based strains to ground motions conversion method has been proposed³⁷, and is modified and adapted here for real-time processing.
The approaches presented in this manuscript build on the above-mentioned advancements in physics-based EEW<sup>24</sup> and DAS earthquake data processing<sup>37</sup>. The potential of the modified strains to ground motions conversion and the new magnitude expression for EEW are examined in conjunction with a theoretical ground motion prediction equation (GMPE)<sup>40</sup>. In the following sections, we present and validate a computationally efficient real-time protocol that relies on straightforward analytical formulations for the analysis of DAS recorded earthquakes. Strain-rates are converted to ground accelerations using a real-time adapted slant-stack approach. Then, earthquake magnitudes are estimated via an analytical expression derived using the Omega-squared source spectra model<sup>41,42</sup>. Subject to the theoretical model, this magnitude expression is applicable to all body wave (P- and S-waves) far-field ground motion recordings. This scheme is applied to several well-coupled fiber segments along different ocean-bottom fibers. Finally, the magnitude is used to predict PGV and PGA away from the hypocenter using a GMPE<sup>24,40</sup>, derived using the same Omega-squared source model. The fact that both magnitude estimation and ground motion prediction are derived from the same theoretical model contributes to the stability and consistency of the estimates, as shown in the next sections. Magnitude and peak ground shaking predictions are continuously updated and modified as new seismic signals are recorded. We demonstrate the robustness of these real-time approaches for a wide magnitude range and show that using well-coupled offshore fibers for EEW can significantly improve warning times compared to those expected from standard point-sensor-based EEW systems.

To demonstrate the merits of the proposed DAS-based EEW schemes, we compiled a DAS earthquake dataset from different tectonic environments. Data were recorded by four different ocean-bottom fibers: two offshore Greece<sup>13,16,37,43</sup>, one offshore France<sup>13,19,37,44</sup> and one offshore Chile (Supplementary Fig. 2). The measurements in Greece, France and Chile were conducted using three different interrogator units: a Febus A1 DAS interrogator, an Aragon Photonics hDAS interrogator and an ASN OptoDAS interrogator, respectively (See “Earthquake dataset” in Methods). We analyzed a total of 53 DAS recorded earthquakes that range from magnitude 2 to 5.7 (Supplementary Fig. 2) at hypocentral distances of 17 to 365 km (Supplementary Fig. 3). Earthquake metadata are provided in Supplementary Table 1.

## Results

**Using DAS data for magnitude estimation.** Ideally, moment magnitude should be estimated using seismic recordings of ground displacements, *D*, rather than ground velocities, *V*, or accelerations, *A*, and the signals should include as much of their low-frequency portion as possible to avoid magnitude saturation<sup>2,3,5</sup>. Ground displacements can be obtained from well-coupled DAS measurements by integrating strain measurements in time (or double integration of strain-rates) and dividing them by the apparent slowness ($n=0$ in Eq. 1)<sup>11,37-39</sup>. This conversion approach has been demonstrated by previous studies that considered DAS instrument response<sup>11</sup> and coupling<sup>45</sup>. However, the use of DAS converted ground displacements is challenging given the inherently high instrumental noise levels, especially at large distances along long fibers<sup>13,46,47</sup>. The behavior of DAS instrumental noise is demonstrated in Fig. (1) for an earthquake of magnitude 3.6 recorded at a distance of 135 km from an optical fiber offshore southeastern Greece<sup>13,16,37,43</sup> (See map in Supplementary Fig. 2). At low frequencies, the instrumental noise of the time-integral of strains ($\propto D$, Fig. 1a,b), strains ($\propto V$, Fig. 1c,d) and strain-rates ($\propto A$, Fig. 1e,f) is proportional to $f^{-2}$, $f^{-1}$ and independent of frequency, $f$, respectively. As a result, strains-integral ($\propto D$, Fig. 1a) and strains ($\propto V$, Fig. 1c) time-series are contaminated by low-frequency noise, and their use may lead to magnitude over-estimation and false alarms. Thus, we only use strain-rates ($\propto A$) for real-time magnitude estimation even though they present a weaker correlation with earthquake magnitude compared to strains-integral ($\propto D$) and strains ($\propto V$) (See “The relation between earthquake source parameters and ground motions” in Methods). Since strain-rates’ instrumental noise increases as $f$ at high frequencies (Fig. 1f), a lowpass filter is needed. This filter will not bias magnitude estimations because larger earthquakes produce lower frequency radiation.

**Strain-rates to ground accelerations conversion.** The performance of the conversion algorithm (See “Real-time strain-rates to ground accelerations conversion” in Methods) is demonstrated for a magnitude 3.8 earthquake recorded offshore Chile at a hypocentral distance of 60 km (See map in Supplementary Fig. 2) for a single DAS channel at a distance of 103 km along the fiber (Fig. 2). Note that direct P-waves are not visible, although P-wave induced scattered-waves are clearly seen (1–6 s in Fig. 2a,b). The same analysis for the largest earthquake in the dataset, a magnitude 5.7, is shown in Supplementary Fig. 4; for this earthquake, strain-rate amplitudes exhibit some saturation (See Discussion). The real-time slant-stack approach resolves the apparent velocities of the different seismic phases: ~4.2 km/s for direct S-waves (6–9 s in Fig. 2a,b) and ~1.8 km/s for surface-waves (e.g., 1–6 s and 10–13 s in Fig. 2a,b). Owing to these velocity variations, ground accelerations are somewhat different from strain-rates: accelerations (blue curve in Fig. 2c) exhibit a noticeable amplitude difference between fast S-waves and slow surface-waves, while strain-rates (black curve in Fig. 2c) display similar amplitudes for both phases. A comparison between the performance of the real-time slant-stack conversion and the previously presented approach<sup>37</sup> indicates that the real-time adaptations do not decrease the conversion quality (Supplementary Fig. 5).

**The effect of stress drop variability.** Stress drop, $\Delta\tau$, is a fundamental earthquake source parameter that strongly affects ground motion intensities<sup>40,48-50</sup> (See “The relation between earthquake source parameters and ground motions” in Methods). For optimal ground motion prediction, both magnitude and stress drop should be determined, as demonstrated by recent studies<sup>3,36,40</sup>. Since in this framework we only use one ground motion metric, i.e., ground accelerations rms, $A_{rms}$, we may only estimate the magnitude (see “Magnitude estimation from bandlimited ground accelerations” in Methods) while the stress drop needs to be set a priori. Because $A_{rms}$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** DAS instrumental noise. (a,b) Strains-integral, (c,d) strains, and (e,f) strain-rates of a magnitude 3.6 earthquake recorded at a hypocentral distance of 135 km at 21 km along a fiber offshore Greece. Signals filtered between 0.06 and 10 Hz are shown on the left as a function of time from P-wave arrival (a,c,e) and the corresponding spectra are shown on the right (blue curves b,d,f). The prefiltered spectra (orange curves b,d,f) demonstrate that low frequency noise is **b** ∝ *f*⁻², **d** ∝ *f*⁻¹, and *f* independent of frequency.</figcaption></figure>

are highly affected by the stress drop, and because its a priori value may deviate from its earthquake-specific
real value⁴⁰,⁵¹, it is useful to examine the effect of stress drop variability on magnitude estimation and intense
shaking prediction. To this end, we synthesized A<sub>rms</sub> using an ideal lowpass Butterworth filter, and PGV<sub>synt</sub> and
PGA<sub>synt</sub> for different magnitudes using Δτ = 10 MPa at a hypocentral distance of 50 km (See “Synthetic ground
motions” in Methods). We then used the synthetic A<sub>rms</sub> to estimate the magnitudes, using different a priori stress
drops of 1, 10 and 100 MPa (Eq. 7). The estimated magnitude and a priori stress drop were then used to predict
PGV<sub>pred</sub> and PGA<sub>pred</sub> (Eq. 10), assuming that the distance is known (Fig. 3). When using Δτ = 10 MPa in Eqs. (7)
and (10), magnitude, PGV, and PGA discrepancies are small (panels b, d and f of Fig. 3, respectively) and mostly
attributed to the approximations made in deriving Eq. (7) (See Supplementary Note 1). When the stress drop in
Eq. (7) is under-estimated (Δτ = 1 MPa) and over-estimated (Δτ = 100 MPa), magnitudes are over-estimated and
under-estimated, respectively, by as much as 1.33 magnitude units for large earthquakes (Fig. 3a). When these
biased magnitudes and stress drops are used to predict PGV, and PGA, they result in reasonable predictions:
the standard deviations of the residuals are limited to ~0.43 log₁₀(PGV) and log₁₀(PGA) units (solid curves in
Fig. 3c,e, respectively). This behavior is explained by inspecting Eq. (10): to first order⁴⁰, PGV ∝ M₀¹/² Δτ¹/² and
PGA ∝ M₀¹/³ Δτ²/³, so using under-estimated stress drops along with over-estimated seismic moments (and
vice-versa), as is the case here, will result in relatively small PGV and PGA discrepancies; Magnitude and stress
drop biases reduce each other's effect on ground motion predictions. In contrast, if synthetic magnitudes are
used in conjunction with the over- and under-estimated stress drops, PGV and PGA discrepancies would be
significantly higher (dashed curves in Fig. 3c,e). Further explanations on the shape of the residual plots are pro-
vided in Supplementary Note 2.

When implementing the proposed methods to different fibers in different tectonic settings, a priori stress
drop may be estimated using available earthquake observations⁴⁰,⁵²,⁵³ or taken from previous studies, if available.
However, the results in Fig. 3 show that while the discrepancies between the synthetic earthquake stress drop
and that used in Eq. (7) may have a significant impact on magnitude estimation, the effect on ground motion
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Strain-rates to ground accelerations conversion and magnitude estimation. **a** Strain-rates of a magnitude 3.8 earthquake recorded at a hypocentral distance of 60 km between 102.3 km and 103.6 km along a fiber offshore Chile. The fiber segment used for magnitude estimation is color-coded (102.7 km to 103.3 km). **b** Semblance as functions of apparent slowness and time from P-wave arrival for a reference DAS channel at 103 km from the interrogator (black dashed line in **a**). Smoothed slowness (See Methods) is indicated by a red curve. **c** Strain-rates (black) and converted ground accelerations (blues) for the reference DAS channel. **d** Real-time magnitude evolution using stress drops of 1 MPa (dashed curve) and 10 MPa (solid curve). Manually picked P- and S-wave arrival times are indicated by magenta lines.</figcaption></figure>

prediction is minimized, and the approach may be reliably used even with a biased stress drop. The effect of stress drop variability will be further examined using recorded earthquakes in the following section.

**Real-time magnitude estimation and peak ground shaking prediction.** The performance of the real-time strain-rates to ground accelerations conversion, magnitude estimation, and ground motion prediction are demonstrated using a composite earthquake catalog of 53 DAS and point-sensor (seismometer and accelerometer) recorded earthquakes from Greece, France, and Chile (See "Earthquake dataset" in Methods, earthquake catalog in Supplementary Table 1, and maps showing the locations of earthquakes, fibers, and point-sensors in Supplementary Fig. 2). These earthquakes range from magnitude 2 to 5.7 (Supplementary Fig. 3) and were recorded by four different offshore fibers using three different DAS interrogators. DAS records are converted to ground accelerations and used to estimate the magnitude, while point-sensor records are used to compare observed and predicted PGV and PGA. Earthquake locations (and hypocentral distances) and P- and S-wave arrival times are assumed to be known: the former are extracted from available earthquake catalogs and
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The effect of stress drop variability on magnitude estimation and ground motion prediction. Estimated minus synthetic magnitude as a function of synthetic magnitude for (a) under-estimated (1 MPa), over-estimated (100 MPa) and (b) known (10 MPa) stress drops. The logarithms of predicted peak ground motions minus the logarithms of synthetic peak ground motions as functions of synthetic magnitude are shown for PGV for (c) 1 MPa and 100 MPa and (d) 10 MPa, and for PGA for (e) 1 MPa and 100 MPa and (f) 10 MPa. The effect of using synthetic magnitude on (c) PGV and (e) PGA discrepancies is indicated by semi-transparent dashed curves. In all panels, curves corresponding to 1, 10 and 100 MPa are indicated by red, black and blue curves, respectively.</figcaption></figure>

the latter are manually picked. In practice, earthquake location and phase picking will be achieved in real-time via additional modules, whose development is beyond the scope of this manuscript. Thus, the uncertainties and discrepancies reported in this section are expected to be larger when earthquake detection and location are also implemented in real-time.

As previously stated, the goal of an EEW system is to produce robust ground motion predictions, while magnitude estimations are merely a by-product. In addition, while we estimate moment magnitudes, most catalogs report local magnitudes, whose values may significantly differ⁵⁴. Thus, in subsequent analysis, we focus on the discrepancies between predicted and observed PGV and PGA as a measure for the algorithms' performance, and provide less attention to the agreement between real-time and catalog reported magnitudes.

Magnitude is estimated using several well-coupled fiber segments for each cable. Coupling quality is evaluated by inspecting earthquakes' seismic wavefield along the fiber and identifying sections that display continuous seismic wavefronts and small amplitude variabilities of less than 4 dB³⁷. DAS data are converted to ground accelerations and an initial magnitude estimate is obtained two seconds following P-wave detection at the first fiber segment, and is continuously updated with increasing data intervals and as the earthquake is recorded at additional locations along the fiber. The analysis uses all available phases including direct P- and S-arrival, scattered waves and surface waves. For P-waves, which are seldom undetected by horizontal DAS arrays, scattered and later arriving P-phases are used for the analysis. In this work, phases were identified and picked manually, while in real-time it will be achieved via automatic algorithms²¹,³¹. Real-time magnitude estimation is demonstrated for a catalog magnitude 3.8 earthquake using a single fiber segment in Fig. 2d. Magnitude increases with time, starting at the scattered P-waves (2–7 s), followed by a significant increase with the S-wave arrivals (7–9.5 s). As theoretically predicted (Fig. 3), magnitude estimates vary for different a priori stress drops, with magnitudes of 5.8 and 4.6 for 1 and 10 MPa, respectively, at 9.5 s from P-wave detection. Similar behavior is observed for the catalog Magnitude 5.7 earthquake shown in Supplementary Fig. 4d. Magnitude estimates improve with time as seen in Fig. 4a–c where real-time and catalog magnitudes are compared at 4, 10 and 15 s from the first P-wave detection, for the entire dataset.

A comparison between predicted (Eq. 10) and observed (See “Earthquake dataset” in Methods) PGV and PGA at 15 s indicates that the residuals are independent of hypocentral distance (Fig. 4d,e and Supplementary Fig. 6d, e) and catalog magnitude (Supplementary Fig. 7d, e), and that their standard deviations are relatively small, only slightly higher than the optimal values, i.e., within-event variabilities reported in the caption of Fig. 4.
<figure><img src="image_6.png" /></figure>

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

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

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

<figure><img src="image_10.png" /><figcaption>**Figure 4.** Real-time magnitude estimation and ground motion prediction using 10 MPa. Real-time magnitude as a function of catalog magnitude at a 4, b 10 and c 15 s from the P-wave arrival at the first fiber segment. Fiber-segment-specific estimates and event averages are indicated by black and red symbols, respectively. The dashed black line is a 1:1 line and the standard deviations of the magnitude residuals are indicated in the bottom right corners for segment specific (black) and event averaged (red) estimates. Discrepancies between the logarithms of predicted and observed peak ground motions are plotted for **d** PGV and **e** PGA as functions of hypocentral distance. Color-code corresponds to catalog magnitudes. Earthquakes from Chile, Greece and France are indicated by stars, circles and triangles, respectively. Panel legends indicate the following: cable name (number of earthquakes, number of PGV and PGA observations), (average residuals, standard deviation to the residuals). Average within event variabilities, i.e., the optimal standard deviation to the residuals, for PGV are 0.68, 0.5 and 0.52 for Chile, Greece and France, respectively, and for PGA are 0.71, 0.61 and 0.59 for Chile, Greece and France, respectively.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 5.** Time-gain using DAS offshore Chile. The region is indicated by a black rectangle in the inset map in **b**. **a** Catalog ($M_{cat}$) and real-time ($M_{RT}$) magnitude estimates when S-waves are expected to reach the coastline. Earthquakes are indicated by circles with size corresponding to catalog magnitudes, and color corresponding to $M_{RT}$ estimation times. The shortest path to the coastline is indicated by grey dashed lines for each earthquake. The fiber is indicated by a blue curve and fiber segments used for magnitude estimation are indicated by black rectangles. **b** P-wave time-gain (red color scale) for different possible earthquake locations. Only positive time-gains are shown; negative time gains indicate earthquake locations that are closer to the fourth closest point-sensors than to the fiber. Point-sensors are indicated by blue triangles. The region shown in **a** is indicated by a black dashed rectangle in **b**. Maps were generated using Python's Basemap package and bathymetric data downloaded from ncei.noaa.gov/maps/bathymetry/.</figcaption></figure>

The latter result suggests that peak ground motion residuals are mainly caused by different site and path conditions that may be accounted for in future implementations, subject to additional research. While magnitude estimates are highly sensitive to the a priori stress drop, PGV and PGA residuals exhibit low sensitivity (Fig. 4 and Supplementary Fig. 6 for 10 and 1 MPa, respectively). This behavior is further demonstrated by examining the average magnitude, and PGV and PGA residuals for the largest available earthquake (Supplementary Fig. 8): average residuals show little sensitivity to stress drop and similar trends to those theoretically predicted (Fig. 3), i.e., PGV residuals are higher for stress drop under-estimation, and PGA residuals are lower for stress drop under-estimation.

## Discussion

The results presented in this manuscript demonstrate that DAS can be reliably used for real-time magnitude estimation and ground motion prediction, two critical components of an operational EEW system. The use of DAS for EEW presents several significant advantages compared to the use of standard point-sensors, especially in the time-gain for offshore earthquakes. This latter advantage is illustrated in Fig. 5 using the fiber deployed offshore Chile, where ocean-bottom earthquakes pose a significant seismic hazard. For the offshore earthquakes shown in Fig. 5a, by the time S-waves are expected to reach the Chilean coastline, real-time magnitude estimates are typically within half a magnitude unit of catalog values, allowing for robust alert issuance before intense ground shaking is felt onland, and well before earthquakes are recorded by the available seismic network (Fig. 5b). The time-gain achieved by using the offshore Chile fiber compared to the current point-sensor network is defined here as the difference between the P-wave arrival at the closest fiber segment and at the fourth seismic station, as commonly required by EEW systems⁵⁵. This time-gain may be as large as 25 s for earthquakes that occur near the fiber and may even result in early detection and alert issuance for onland earthquakes where point-sensor
coverage is sparse (Fig. 5b). These precious seconds can have a decisive impact on mitigating the risk posed by potentially catastrophic offshore earthquakes.

Together with the time-gain for offshore earthquakes, DAS-based EEW can potentially outperform point-sensor-based EEW for several reasons. When implemented on well-coupled fiber segments, magnitude estimates are more reliable since data from many closely spaced DAS channels are averaged, reducing the impact of outliers and smoothing local effects. DAS facilitates robust differentiation between earthquakes and noise since earthquakes' seismic wavefield is near-instantaneously recorded on hundreds-of-meters long fiber segments. As a result, false detections will be reduced and one fiber segment may be sufficient to issue early warning, subject to earthquake location capabilities.

While direct S-waves are detected by horizontally installed fibers, direct P-waves are usually not (Fig. 2a), a result of their fast velocities and the angle between the waves' polarization and the fiber's axis<sup>13,56</sup>. In contrast, P-wave induced scattered-waves are well recorded (2–6 s in Fig. 2a,b). Direct P-waves, if available, and S-waves, as well as scattered P- and S-waves are all used for magnitude estimation. While scattering results from heterogeneities of Earth's media and varies from one region to another, the use of the waves' apparent velocities to convert strain-rates to ground accelerations reduces the effect of this local phenomena on magnitude estimations. The dominance of these scattered waves will pose difficulties for earthquake location, since using scattered P-waves instead of direct P-waves will likely point to the scatterers' locations rather than the earthquake's source. Because P-wave based magnitudes are typically under-estimated (Fig. 2d), they are not expected to cause false alarms, yet they may be sufficient to surpass predefined alert thresholds. Since for EEW, sensors should be placed at proximity to expected epicenters, the closer the fiber is to earthquake locations, the sooner the high signal-to-noise S-waves are detected and used. For large earthquakes, the lower sensitivity of DAS to P-waves is an advantage because it will limit the saturation of direct and scattered P-waves.

Since DAS is an emerging technology, available datasets do not contain sufficient earthquakes whose damage potential is of interest for EEW. As a result, several technical aspects of the proposed schemes such as amplitude saturation and fiber-ground coupling during strong shaking cannot be fully evaluated. However, unlike commonly used empirical EEW approaches<sup>32–34</sup>, the proposed scheme is theoretical and relies on a well-established source model<sup>41</sup> that was found to adequately describe the far-field radiation of a large range of earthquakes. Thus, showing that the proposed methods work for the current earthquake dataset is sufficient to demonstrate their validity. In addition, considering small magnitude earthquakes, as we do here, demonstrates the robustness of the system to false alarms. Analyzing near field records is a troublesome issue for EEW; since a complete theoretical framework is yet to be developed<sup>57</sup>, other approaches such as resolving line sources<sup>58</sup> or extrapolating peak ground motions away from the earthquake source<sup>59</sup> may need to be adapted to DAS data in order to address this gap.

The derivation of the presented physics-based magnitude estimation approach did not require any earthquake observations, a significant advantage since the scarcity of DAS earthquake observations hinders the derivation of empirical methods. Because no earthquakes were required, the approach is geographically independent and readily applicable in any tectonic setting using both offshore and onland fibers and different DAS interrogator units, as demonstrated here using earthquakes from Greece, France and Chile. Earthquake observations are only required to map well-coupled fiber segments, although this objective may also be achieved using ambient noise<sup>16</sup>. Using segments that are not well-coupled may lead to either magnitude under-estimation, if strain amplitudes are weak, or over-estimation, if a segment is suspended and experiences strong vibrations due to cable-waves<sup>60</sup>. The approach allows for continuous update of magnitude and ground motion predictions, key for analyzing large earthquakes with long durations. In addition, using a holistic magnitude estimation and ground motion prediction that are derived from the same earthquake model reduces the impact of stress drop related magnitude biases on ground motion predictions and enhances the overall robustness of the system.

The computational costs of the presented approaches are low. While DAS acquisitions typically provide very large data volumes, for EEW, data can be largely down-sampled in both time and space, limiting both the volume of data and processing times. For instance, to obtain timely and robust magnitude estimates, it is sufficient to use preselected well-coupled fiber portions at spacings of several kilometers. For the purpose of this study, we analyzed 180 s long recordings of preselected DAS fiber segments (33 channels), down-sampled to ~20 Hz, in ~136 s using a Python code on an Intel Core i7 laptop with 32 GB RAM using a single thread. These computation times indicate that the method is valid for real-time. For future implementation as part of an operational EEW system, several aspects of the code can be optimized and run in parallel. For example, slant-stack computations, which are the most time-consuming (~4 s per channel), can be parallelized, in addition to computations for different fiber segments. The latter can also increase the number of fiber segments used for magnitude estimation.

For few earthquakes, strain amplitudes exhibited a small degree of saturation. Nevertheless, magnitude estimations still allow for reliable ground motion predictions (Fig. 4). This phenomenon is insufficiently reported and investigated in existing literature and needs to be quantified and addressed as it may affect the ability to analyze higher strain amplitudes and provide reliable warnings for larger earthquakes. DAS saturation needs to be studied along with DAS manufacturers to devise preprocessing and postprocessing methods in order to fully demonstrate the viability of DAS for EEW.

The framework presented in this study demonstrates the great potential of using DAS for EEW. The approaches presented here allow for easy, robust, and fast implementation of EEW using both offshore and onland optical fibers in any tectonic setting. Specifically, using existing ocean-bottom optical fibers, which are almost ubiquitous along subduction zones worldwide, provide a cheap and readily available EEW solution, especially for exposed developing countries, that will significantly enhance earthquake hazard mitigation capabilities.
Methods

**The relation between earthquake source parameters and ground motions.** For large earthquakes, i.e., when high-frequency attenuation is negligible, recorded in the far-field, ground displacements root-mean-squares (rms), D<sub>rms</sub>, and peak ground displacements (PGD) are mostly a function of the seismic moment, M<sub>0</sub>:

D<sub>rms</sub> ~ PGD ~ M<sub>0</sub><sup>5/6</sup> Δτ<sup>1/6</sup> while ground velocities rms, V<sub>rms</sub>, and PGV, and accelerations rms, A<sub>rms</sub>, and PGA, are also strongly influenced by the stress drop, Δτ: V<sub>rms</sub> ~ PGV ~ M<sub>0</sub><sup>1/2</sup> Δτ<sup>1/2</sup> and A<sub>rms</sub> ~ PGA ~ M<sub>0</sub><sup>1/3</sup> Δτ<sup>2/3</sup> <sup>24,40</sup>. The proportionality between rms and peak ground motions stems from statistical theories<sup>61</sup> and was observed by previous studies<sup>40</sup>. Note the different powers associated with M<sub>0</sub> and Δτ. Thus, ground displacements serve as a better magnitude predictor compared to velocities or accelerations<sup>32,35,40</sup>.

**Real-time strain-rates to ground accelerations conversion.** The slant-stack<sup>62</sup> based strains to ground motions conversion scheme<sup>37</sup> accounts for apparent phase velocity variations in both time and space. The conversion is applied for each DAS channel along the fiber using short, approximately linear, fiber segments. Here, this recently presented approach<sup>37</sup> is modified and optimized for real-time performance. The semblance (coherency) as a function of apparent slowness p<sub>x</sub> and time t, for a DAS channel located at x<sub>0</sub> along the fiber, can be written as:

$$sem(p_x, t) = \frac{1}{L} \begin{cases} \frac{\left[ \sum_{j=-L}^{L} g(t+p_x(x_j-x_0)) \right]^2}{\sum_{j=-L}^{L} g(t+p_x(x_j-x_0))^2} & \text{if } p_x > 0 \\ \frac{\left[ \sum_{j=1}^{L} g(t+p_x(x_j-x_0)) \right]^2}{\sum_{j=1}^{L} g(t+p_x(x_j-x_0))^2} & \text{if } p_x < 0 \end{cases}, \quad (2)$$

where L is the number of DAS channels used for slowness estimation, g(t) is the DAS strain-rates time-series, and x<sub>j</sub>—x<sub>0</sub> is the distance between station j and the reference channel (at x<sub>0</sub>). Equation (2) can be regarded as the causal slant-stack, where only data samples of g(t) that have already been recorded are considered.

The conversion procedure is performed as follows. For computational efficiency, recorded strain-rates are down-sampled to 20 Hz (or slightly higher, depending on the original signals' sampling-rate). Data is lowpass filtered at 5 Hz using a 4-pole Butterworth filter to diminish high frequency instrumental noise. The applied down-sampling and filtering did not decrease the robustness of the conversion and subsequent magnitude estimation. The local slant-stack transform is applied using fiber segments of ~380 m length<sup>37</sup>, with channel spacings of ~20 m, skipping several channels for densely spaced measurements. The used fiber segments are long enough to resolve long seismic wavelengths with fast velocities of several km/s, and short enough so that seismic waves are coherent and fiber sections are approximately linear<sup>37</sup>. Semblance is calculated using 50 predefined slowness values, equally spaced between -5 s/km and 5 s/km. At each t, the wavefield's slowness is determined as the one with highest semblance. The produced slowness time-series is then smoothed by applying a causal moving-mean filter of 1 s to its absolute value. Strain-rates time-series are then divided by the slowness time-series to obtain ground accelerations, followed by an additional 5 Hz lowpass filter. Because we are eventually interested in the converted strain-rates' rms, the slowness' sign may be discarded (See "Magnitude estimation from bandlimited ground accelerations" in Methods).

**Magnitude estimation from bandlimited ground accelerations.** We derive an expression for the rms of the ground accelerations using the commonly used Omega-squared source model<sup>41</sup> describing the far-field body wave spectra (grey dashed curve in Supplementary Fig. 9). This derivation procedure follows that used by several previous studies<sup>24,35,40,63–66</sup>. The acceleration omega-squared model<sup>41</sup> subject to high frequency attenuation<sup>67</sup> (grey dotted curve in Supplementary Fig. 9) reads as:

$$\ddot{\Omega}(f) = (2\pi f)^2 \frac{\Omega_0}{1 + \left(\frac{f}{f_0}\right)^2} e^{-\pi \kappa f}, \quad (3)$$

where f<sub>0</sub> is the source corner frequency, Ω<sub>0</sub> is the displacement low frequency spectral plateau, and k is an attenuation parameter. Since strain-rates are lowpass filtered at 5 Hz, acceleration rms are calculated using Eq. (3) as

A<sub>rms</sub> = √(2/T) ∫<sub>f=0</sub><sup>f=5</sup> |Ω̄(f)|<sup>2</sup> df (black dashed curve in Supplementary Fig. 9), where T is the data interval. The integral is solved and an analytic approximation is obtained (See Supplementary Note 1). The spectral parameters Ω<sub>0</sub> and f<sub>0</sub> are substituted with the seismic moment<sup>23</sup> and stress drop<sup>68</sup>, respectively, via:

$$M_0 = \Omega_0 \frac{4\pi \rho C^3 R}{U_{\phi\theta} F_s}, \quad (4a)$$

$$\Delta\tau = \frac{7}{16} M_0 \left( \frac{f_0}{kC_s} \right)^3, \quad (4b)$$

where ρ is the density at the source, C is the wave velocity at the source (C<sub>p</sub> and C<sub>s</sub> for P- and S-waves, respectively), R is the hypocentral distance, U<sub>φθ</sub> is the average radiation pattern, F<sub>s</sub> is the free-surface correction, and k is a phase-specific constant which also depends on the source model and rupture speed<sup>42</sup>. Equation (4b) is valid for a circular crack embedded in a homogeneous medium<sup>68</sup>. The resulting expression is:
$$A_{rms}^{approx} = A_1 M_0^{\frac{1}{3}} \Delta \tau^{\frac{2}{3}} \sqrt{1 - e^{-2\alpha_m}} \frac{1}{R\sqrt{\kappa T} \left( 1 + \frac{A_2^2 \kappa^2 (\frac{\Delta \tau}{M_0})^{\frac{2}{3}} \sqrt{1 - e^{-2\alpha_m}}}{h(\alpha_m)} \right)}, \quad (5)$$

where the superscript *approx* signifies approximate rms, $A_1 = \frac{U_{\varphi\theta} F_s \sqrt{\pi}}{\rho C^3} (\frac{16}{7})^{\frac{2}{3}} (kC_s)^2$, $A_2 = \pi (\frac{16}{7})^{1/3} kC_s$, $h(\alpha_m) = e^{-\alpha_m} \sqrt{\frac{1}{2}[-3 - 6\alpha_m - 6\alpha_m^2 - 4\alpha_m^3 - 2\alpha_m^4 + 3e^{2\alpha_m}]}$ and $\alpha_m = 5\pi\kappa$.

Equation (5) can be analytically solved for the seismic moment:

$$M_0 = \frac{1}{27a_1} \left( \frac{a_4}{2^{\frac{1}{3}}} + \frac{2^{\frac{1}{3}}a_2^2}{a_4} + a_2 \right)^3, \quad (6)$$

where $a_1 = A_1 \Delta \tau^{\frac{2}{3}} \sqrt{1 - e^{-2\alpha_m}} \frac{1}{R\sqrt{\kappa T}}$, $a_2 = A_{rms}$, $a_3 = \frac{A_{rms} A_2^2 \Delta \tau^{\frac{2}{3}} \kappa^2 \sqrt{1 - e^{-2\alpha_m}}}{h(\alpha_m)}$ and $a_4 = (3\sqrt{3}(27a_1^4 a_3^2 + 4a_1^2 a_2^3 a_3) + 27a_1^2 a_3 + 2a_2^3)^{\frac{1}{3}}$. The moment magnitude can then be written as:

$$M_W = 2 \left( \frac{a_4}{2^{\frac{1}{3}}} + \frac{2^{\frac{1}{3}} a_2^2}{a_4} + a_2 \right) - \frac{2}{3}(a_1) - 7.05, \quad (7)$$

where $M_0$ is expressed in Nm.

While the coefficients $a_1, a_2, a_3$ and $a_4$ contain many parameters, only few are updated in real-time: $A_{rms}$ is continuously updated as new data is recorded, the available data interval $T$ begins at the P-wave arrival and increases with time, and $R$ is updated as earthquake location improves. The parameters used are$^{24}$: $F_s = 2$, $\rho = 2600$ kg/m$^3$, $C_s = 3.2$ km/s, $C_p = 5.3$ km/s, $\kappa = 0.025$ s, $U_{\varphi\theta}$ equals 0.52 and 0.63 for P- and S-waves, respectively$^{23}$, and $k$ equals 0.32 and 0.21 for P- and S-waves, respectively$^{42}$. For data intervals that contain both P- and S-waves, the phase specific constants need to be averaged based on the relative intervals of each phase$^{24}$:

$$\text{const} = \frac{T_{S-P}}{T} \text{const}_P + \frac{T - T_{S-P}}{T} \text{const}_S, \quad (8)$$

where const stands for $U_{\varphi\theta}$, C or k for P- or S-waves, and $T_{S-P}$ is the S-P data interval. Using these parameters, $a_1$ and $a_3$ may be written as:

$$a_1 = 113014 \left( \frac{k^2 U_{\varphi\theta}}{C^3} \right) \Delta \tau^{\frac{2}{3}} \frac{1}{R\sqrt{T}}, \quad (9a)$$

$$a_3 = 1828968(k^2) \Delta \tau^{\frac{2}{3}} A_{rms}, \quad (9b)$$

where phase-specific terms are written in parentheses.

In this application, the magnitude is estimated using several manually identified well-coupled fiber segments of ~600 m as follows. Strain-rates within each fiber segment are converted to ground accelerations (See “Real-time strain-rates to ground accelerations conversion” in Methods). $A_{rms}$ is calculated per DAS channel starting at the P-wave arrival, and is then logarithmically averaged per fiber segment at every time-step to minimize the impact of outliers. Since DAS can only measure the wavefield in-line with the fiber, $A_{rms}$ is multiplied by $\sqrt{2}$ to compensate for the missing orthogonal component. The averaged $A_{rms}$ at time $T$ is then input to Eq. (7) along with $\Delta\tau$ and $R$ to estimate the magnitude. Magnitude estimates are continuously updated until either averaged $A_{rms}$ reaches its maximum value, or $T=60$ seconds$^{24}$. Magnitude estimates from different fiber segments are weight-averaged by the available data interval to obtain an event specific estimate.

**Ground motion prediction.** For PGV and PGA prediction, we use a set of physics-based GMPEs$^{24,40}$, derived using the same source model$^{41}$ (Eq. 3) used to obtain the real-time magnitude expression (Eq. 7) (See “Magnitude estimation from bandlimited ground accelerations” in Methods). The GMPEs for PGV and PGA are:

$$\text{PGV} = 2.9 \sqrt{M_0 \Delta \tau} \frac{\beta_V}{R \sqrt{\frac{1}{kC_s} \left( \frac{7}{16} \frac{M_0}{\Delta \tau} \right)^{1/3} + R/C_s} \left[ 1 + \pi^{\frac{4}{3}} \kappa_0 k C_s \left( \frac{16}{7} \frac{\Delta \tau}{M_0} \right)^{\frac{1}{3}} \right]^{\frac{3}{2}}} \quad (10a)$$

$$\text{PGA} = 3.3 M_0^{1/3} \Delta \tau^{2/3} \frac{\beta_A}{R \sqrt{\kappa_0 \left[ \frac{1}{kC_s} \left( \frac{7}{16} \frac{M_0}{\Delta \tau} \right)^{\frac{1}{3}} + \frac{R}{C_s} \right] \left[ 1 + 1.5^{-\frac{1}{4}} \pi \kappa_0 k C_s \left( \frac{16}{7} \frac{\Delta \tau}{M_0} \right)^{\frac{1}{3}} \right]^2}}, \quad (10b)$$
where $\beta_V = \frac{2\pi U_{\phi\theta} F_s \sqrt{\frac{16}{7} (kC_S)^{\frac{3}{2}}}}{(\sqrt{2\pi 4\rho C_S^3})}$ and $\beta_A = \frac{4\pi U_{\phi\theta} F_s (\frac{16}{7})^{2/3} (kC_S)^2}{(\sqrt{\pi 4\rho C_S^3})}$. These theoretical GMPEs are readily applicable in any seismic region. Using the parameter tuning for S-waves (See "Magnitude estimation from bandlimited ground accelerations" in Methods), $\beta_V = 2.44 \times 10^{-10} \text{ m}^{1.5}\text{s}^{1.5}/\text{kg}$ and $\beta_A = 2.05 \times 10^{-8} \text{ m}^2\text{s}/\text{kg}$.

**Synthetic ground motions.** The GMPEs in Eq. (10) are used to generate synthetic PGV and PGA for different seismic moments, stress drops and hypocentral distances. Synthetic $A_{rms}$ are generated by calculating the rms of the acceleration spectra (Supplementary Fig. 9). These spectra are produced for a specific seismic moment, stress drop and hypocentral distance using Eq. (3) and (4), subject to a lowpass filter. The filter is modeled in two manners: as a clean cutoff (dashed black curve in Supplementary Fig. 9) as used for the model derivation (See "Magnitude estimation from bandlimited ground accelerations" in Methods), or as an ideal 4-pole Butterworth filter (solid black curve in Supplementary Fig. 9), similar to that used for DAS signal processing.

**Earthquake dataset.** DAS measurements in Greece were conducted using a Febus A1 DAS interrogator between 18–19 and 19–25 April 2019 on 13.2 km and 26.2 km long fibers, sampled at 6 ms and 5 ms, respectively. Gauge length and spatial sampling were both set to 19.2 m for the two fibers. DAS measurements in France were conducted using an Aragon Photonics hDAS interrogator between 11–31 July 2019 on a 44.8 km long fiber, sampled at 10 ms and 2 ms for the first and last 10 days, respectively. Gauge length and spatial sampling were both set to 10 m. DAS measurements in Chile were conducted using an ASN OptoDAS interrogator between 27 October and 3 December 2021 on a 204 km long fiber, sampled at 8 ms. Gauge length and spatial sampling were both set to 4.085 m. The Febus and OptoDAS interrogators record strain-rates while the Aragon instrument records stains; the latter were differentiated to strain-rates before the conversion to ground accelerations.

Seismometer and accelerometer recordings were used to calculate PGV and PGA for the different earthquakes as follows. Data for Greece, France and Chile were obtained from the National Observatory of Athens, the RESIF repository and IRIS, respectively. The two horizontal components of point-sensors were demeaned and highpass filtered at 1 Hz using a 4-pole Butterworth filter, followed by a simple gain correction. Velocity-meter signals were differentiated to obtain ground accelerations and accelerometer records were integrated to obtain ground velocities. An additional highpass filter was applied after differentiations and integrations. PGV (PGA) were then calculated as the geometric mean of the maximum of the absolute value of the two velocity (acceleration) components. PGV and PGA that are smaller than 5 times the standard deviation of the associated time-series are discarded as they may be biased by noise.

**Data availability**

Samples of DAS earthquakes are available on https://osf.io/4bjph/.

Received: 23 September 2022; Accepted: 2 January 2023

Published online: 09 January 2023

**References**

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. Colombelli, S., Carotenuto, F., Elia, L. & Zollo, A. Design and implementation of a mobile device app for network-based earthquake early warning systems (EEWSs): application to the PRESTo EEWS in southern Italy. *Nat. Hazards Earth Syst. Sci.* **20**, 921–931 (2020).

3. Nof, R. N., Lior, I. & Kurzon, I. Earthquake early warning system in Israel—Towards an operational stage. *Front. Earth Sci.* **9** (2021).

4. Licciardi, A., Bletery, Q., Rouet-Leduc, B., Ampuero, J.-P. & Juhel, K. Instantaneous tracking of earthquake growth with elastogmetry signals. *Nature* https://doi.org/10.1038/s41586-022-04672-7 (2022).

5. Kuyuk, H. S. & Allen, R. M. Optimal seismic network density for earthquake early warning: A case study from California. *Seismol. Res. Lett.* **84**, 946–954 (2013).

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

7. Aoi, S. et al. MOWLAS: NIED observation network for earthquake, tsunami and volcano. *Earth Planets Space* **72**, 126 (2020).

8. Schlesinger, A. et al. An Earthquake early warning system for Southwestern British Columbia. *Front. Earth Sci.* **9** (2021).

9. Zhan, Z. Distributed acoustic sensing turns fiber-optic cables into sensitive seismic antennas. *Seismol. Res. Lett.* **91**, 1–15 (2020).

10. Fernandez-Ruiz, M. R. et al. Seismic monitoring with distributed acoustic sensing from the near-surface to the deep oceans. *J. Light. Technol.* **40**, 1453–1463 (2022).

11. Lindsey, N. J., Rademacher, H. & Ajo-Franklin, J. B. On the broadband instrument response of fiber-optic DAS arrays. *J. Geophys. Res. Solid Earth* **125**, 1251–1263 (2020).

12. van den Ende, M. P. A. & Ampuero, J.-P. Evaluating seismic beamforming capabilities of distributed acoustic sensing arrays. *Solid Earth* **12**, 915–934 (2021).

13. Lior, I. et al. On the detection capabilities of underwater distributed acoustic sensing. *J. Geophys. Res. Solid Earth* **126** (2021).

14. Yuan, S., Lellouch, A., Clapp, R. G. & Biondi, B. Near-surface characterization using a roadside distributed acoustic sensing array. *Lead. Edge* **39**, 646–653 (2020).

15. Williams, E. F. et al. Scholte wave inversion and passive source imaging with ocean-bottom DAS. *Lead. Edge* **40**, 576–583 (2021).

16. Lior, I., Mercerat, E. D., Rivet, D., Sladen, A. & Ampuero, J. Imaging an underwater basin and its resonance modes using optical fiber distributed acoustic sensing. *Seismol. Res. Lett.* **93**, 1573–1584 (2022).

17. Spica, Z. J. et al. Subsurface imaging with ocean-bottom distributed acoustic sensing and water phases reverberations. *Geophys. Res. Lett.* **49**, 1–8 (2019).

18. Papp, B., Donno, D., Martin, J. E. & Hartog, A. H. A study of the geophysical response of distributed fibre optic acoustic sensors through laboratory-scale experiments: Geophysical response of fibre optic sensors. *Geophys. Prospect.* **65**, 1186–1204 (2017).

19. Sladen, A. et al. Distributed sensing of earthquakes and ocean-solid Earth interactions on seafloor telecom cables. *Nat. Commun.* **10**, 1–8 (2019).

20. Lellouch, A. & Biondi, B. L. Seismic applications of downhole DAS. *Sensors* **21**, 2897 (2021).
21. Li, Z. & Zhan, Z. Pushing the limit of earthquake detection with distributed acoustic sensing and template matching: a case study at the Brady geothermal field. *Geophys. J. Int.* **215**, 1583–1593 (2018).

22. Singh, S., Capdeville, Y. & Igel, H. Correcting wavefield gradients for the effects of local small-scale heterogeneities. *Geophys. J. Int.* **220**, 996–1011 (2020).

23. Aki, K. & Richards, P. G. *Quantitative Seismology* (University Science Books, 2002).

24. Lior, I. & Ziv, A. Generic source parameter determination and ground-motion prediction for earthquake early warning generic source parameter determination and ground-motion prediction for earthquake early warning. *Bull. Seismol. Soc. Am.* **110**, 345–356 (2020).

25. Worden, C. B., Gerstenberger, M. C., Rhoades, D. A. & Wald, D. J. Probabilistic relationships between ground-motion parameters and modified Mercalli intensity in California. *Bull. Seismol. Soc. Am.* **102**, 204–221 (2012).

26. Satriano, C., Lomax, A. & Zollo, A. Real-time evolutionary earthquake location for seismic early warning. *Bull. Seismol. Soc. Am.* **98**, 1482–1494 (2008).

27. Lomax, A., Satriano, C. & Vassallo, M. Automatic picker developments and optimization: FilterPicker-a Robust, broadband picker for real-time seismic monitoring and earthquake early warning. *Seismol. Res. Lett.* **83**, 531–540 (2012).

28. Eisermann, A. S., Ziv, A. & Wust-Bloch, H. G. Array-based earthquake location for regional earthquake early warning: Case studies from the dead sea transform. *Bull. Seismol. Soc. Am.* **108**, 2046–2053 (2018).

29. Li, Z., Meier, M.-A., Hauksson, E., Zhan, Z. & Andrews, J. Machine learning seismic wave discrimination: Application to earthquake early warning. *Geophys. Res. Lett.* **45**, 4773–4779 (2018).

30. Lellouch, A., Yuan, S., Spica, Z., Biondi, B. & Ellsworth, W. L. Seismic velocity estimation using passive downhole distributed acoustic sensing records: examples from the san andreas fault observatory at depth. *J. Geophys. Res. Solid Earth* **124**, 6931–6948 (2019).

31. Nayak, A., Ajo-Franklin, J., & the Imperial Valley Dark Fiber Team. Distributed acoustic sensing using dark fiber for array detection of regional earthquakes. *Seismol. Res. Lett.* **92**, 2441–2452 (2021).

32. Kuyuk, H. S. & Allen, R. M. A global approach to provide magnitude estimates for earthquake early warning alerts: Magnitude for earthquake early warning. *Geophys. Res. Lett.* **40**, 6329–6333 (2013).

33. Sadeh, M., Ziv, A. & Wust-Bloch, H. Real-time magnitude proxies for earthquake early warning in Israel. *Geophys. J. Int.* **196**, 939–950 (2014).

34. Nazeri, S., Shomali, Z. H., Colombelli, S., Elia, L. & Zollo, A. Magnitude estimation based on integrated amplitude and frequency content of the initial P wave in earthquake early warning applied to Tehran, Iran. *Bull. Seismol. Soc. Am.* **107**, 1432–1438 (2017).

35. Lior, I., Ziv, A. & Madariaga, R. P-wave attenuation with implications for earthquake early warning. *Bull. Seismol. Soc. Am.* **106**, 13–22 (2016).

36. Ziv, A. & Lior, I. Real-time moment magnitude and stress drop with implications for real-time shaking prediction. *Bull. Seismol. Soc. Am.* **106**, 2459–2468 (2016).

37. Lior, I. et al. Strain to ground motion conversion of distributed acoustic sensing data for earthquake magnitude and stress drop determination. *Solid Earth* **12**, 1421–1442 (2021).

38. Daley, T. M., Miller, D. E., Dodds, K., Cook, P. & Freifeld, B. M. Field testing of modular borehole monitoring with simultaneous distributed acoustic sensing and geophone vertical seismic profiles at Citronelle, Alabama: Field testing of MBM. *Geophys. Prospect.* **64**, 1318–1334 (2016).

39. Wang, H. F. et al. Ground motion response to an ML 4.3 earthquake using co-located distributed acoustic sensing and seismometer arrays. *Geophys. J. Int.* **213**, 2020–2036 (2018).

40. Lior, I. & Ziv, A. The relation between ground motion, earthquake source parameters, and attenuation: Implications for source parameter inversion and ground motion prediction equations. *J. Geophys. Res. Solid Earth* **123**, 5886–5901 (2018).

41. Brune, J. N. Tectonic stress and the spectra of seismic shear waves from earthquakes. *J. Geophys. Res.* **75**, 4997–5009 (1970).

42. Madariaga, R. Dynamics of an expanding circular fault. *Bull. Seismol. Soc. Am.* **66**, 639–666 (1976).

43. van den Ende, M. et al. A self-supervised deep learning approach for blind denoising and waveform coherence enhancement in distributed acoustic sensing data. *IEEE Trans. Neural Netw. Learn. Syst.* **1**–14 (2021) https://doi.org/10.1109/TNNLS.2021.31328 32.

44. Guerin, G. et al. High resolution ambient noise tomography of the Southwestern Alps and the Ligurian margin. *Geophys. J. Int.* **220**, 806–820 (2020).

45. Ajo-Franklin, J. B. et al. Distributed acoustic sensing using dark fiber for near-surface characterization and broadband seismic event detection. *Sci. Rep.* **9** (2019).

46. Atterholt, J., Zhan, Z., Shen, Z. & Li, Z. A unified wavefield-partitioning approach for distributed acoustic sensing. *Geophys. J. Int.* **228**, 1410–1418 (2021).

47. Fernández-Ruiz, M. R., Costa, L. & Martins, H. F. Distributed acoustic sensing using chirped-pulse phase-sensitive OTDR technology. *Sensors* **19**, 4368 (2019).

48. Oth, A., Miyake, H. & Bindi, D. On the relation of earthquake stress drop and ground motion variability. *J. Geophys. Res. Solid Earth* **122**, 5474–5492 (2017).

49. Trugman, D. T. & Shearer, P. M. Strong correlation between stress drop and peak ground acceleration for recent M 1–4 earthquakes in the San Francisco Bay Area. *Bull. Seismol. Soc. Am.* **108**, 929–945 (2018).

50. Trugman, D. T., Page, M. T., Minson, S. E. & Cochran, E. S. Peak ground displacement saturates exactly when expected: Implications for earthquake early warning. *J. Geophys. Res. Solid Earth* **124**, 4642–4653 (2019).

51. Kaneko, Y. & Shearer, P. M. Variability of seismic source spectra, estimated stress drop, and radiated energy, derived from cohesive-zone models of symmetrical and asymmetrical circular and elliptical ruptures. *J. Geophys. Res. Solid Earth* **120**, 1053–1079 (2015).

52. Trugman, D. T. & Shearer, P. M. Application of an improved spectral decomposition method to examine earthquake source scaling in Southern California: Earthquake source scaling in Southern California. *J. Geophys. Res. Solid Earth* **122**, 2890–2910 (2017).

53. Shearer, P. M., Abercrombie, R. E., Trugman, D. T. & Wang, W. Comparing EGF methods for estimating corner frequency and stress drop from P wave spectra. *J. Geophys. Res. Solid Earth* **124**, 3966–3986 (2019).

54. Deichmann, N. Local magnitude, a moment revisited. *Bull. Seismol. Soc. Am.* **96**, 1267–1277 (2006).

55. Kohler, M. D. et al. Earthquake early warning ShakeAlert 2.0: Pubic Rollout. *Seismol. Res. Lett.* **91**, 1763–1775 (2020).

56. Mateeva, A. et al. Distributed acoustic sensing for reservoir monitoring with vertical seismic profiling: Distributed acoustic sensing (DAS) for reservoir monitoring with VSP. *Geophys. Prospect.* **62**, 679–692 (2014).

57. Madariaga, R., Ruiz, S., Rivera, E., Leyton, F. & Baez, J. C. Near-field spectra of large earthquakes. *Pure Appl. Geophys.* **176**, 983–1001 (2019).

58. Böse, M., Heaton, T. H. & Hauksson, E. Real-time Finite Fault Rupture Detector (FinDer) for large earthquakes. *Geophys. J. Int.* **191**, 803–812 (2012).

59. Kodera, Y. et al. The Propagation of Local Undamped Motion (PLUM) Method: A Simple and Robust Seismic Wavefield Estimation Approach for Earthquake Early Warning. *Bull. Seismol. Soc. Am.* **108**, 983–1003 (2018).

60. Munn, J. D., Coleman, T. I., Parker, B. L., Mondanos, M. J. & Chalari, A. Novel cable coupling technique for improved shallow distributed acoustic sensor VSPs. *J. Appl. Geophys.* **138**, 72–79 (2017).

61. Van Houtte, C., Larkin, T. & Holden, C. On durations, peak factors, and Nonstationarity corrections in seismic hazard applications of random vibration theory. *Bull. Seismol. Soc. Am.* **108**, 418–436 (2018).
62. Taner, M. T., Koehler, F. & Sheriff, R. E. Complex seismic trace analysis. *Geophysics* **44**, 1041–1063 (1979).

63. Hanks, T. C. *b* values and ω⁻¹γ seismic source models: Implications for tectonic stress variations along active crustal fault zones and the estimation of high-frequency strong ground motion. *J. Geophys. Res.* **84**, 2235 (1979).

64. McGuire, R. K. & Hanks, T. C. RMS accelerations and spectral amplitudes of strong ground motion during the San Fernando, California earthquake. *Bull. Seismol. Soc. Am.* **70**, 1907–1919 (1980).

65. Luco, J. E. On strong ground motion estimates based on models of the radiated spectrum. *Bull. Seismol. Soc. Am.* **75**, 641–649 (1985).

66. Lior, I. & Ziv, A. The relation between ground acceleration and earthquake source parameters: Theory and observations. *Bull. Seismol. Soc. Am.* **107**, 1012–1018 (2017).

67. Anderson, J. G. & Hough, S. E. A model for the shape of the fourier amplitude spectrum of acceleration at high frequencies. *Bull. Seismol. Soc. Am.* **74**, 1969–1993 (1984).

68. Eshelby, J. D. The determination of the elastic field of an ellipsoidal inclusion, and related problems. *Proc. R. Soc. Lond. Ser. Math. Phys. Sci.* **241**, 376–396 (1957).

## Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 101041092), from the UCAJEDI Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01, and from the Observatoire de la Côte d'Azur. We thank GTD Grupo SA who provided access to the infrastructure, and el Centro Sismologico Nacional staff who helped in the logistics.

## Author contributions

I.L. designed the presented algorithms, performed the analysis and wrote the initial draft. D.R., J.P.A. and A.S. contributed to the discussion, methodology, interpretation and presentation of the results. D.R., J.P.A., A.S., S.B., R.S.O., G.A.V.O. and J.A.B.P. took part in performing DAS measurements.

## 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/s41598-023-27444-3.

**Correspondence** and requests for materials should be addressed to I.L.

**Reprints and permissions information** is available at 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_12.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 licence, and indicate if changes were made. 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/4.0/.

© The Author(s) 2023