# An analysis of the dynamic range of Distributed Acoustic Sensing for Earthquake Early Warning

Martijn P.A. van den Ende, Alister Trabattoni, Marie Baillet, Diane Rivet

► To cite this version:

Martijn P.A. van den Ende, Alister Trabattoni, Marie Baillet, Diane Rivet. An analysis of the dynamic range of Distributed Acoustic Sensing for Earthquake Early Warning. Seismica, 2025, 4 (1), 10.26443/seismica.v4i1.1371. hal-05171564

HAL Id: hal-05171564

https://hal.science/hal-05171564v1

Submitted on 20 Jul 2025

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>
# An analysis of the dynamic range of Distributed Acoustic Sensing for Earthquake Early Warning

Martijn van den Ende ID * 1, Alister Trabattoni ID 1, Marie Baillet ID 1, Diane Rivet ID 1

¹Université Côte d'Azur, CNRS, IRD, Observatoire de la Côte d'Azur, Géoazur, France

Author contributions: Formal Analysis: Martijn van den Ende. Resources: Marie Baillet, Alister Trabattoni. Writing - original draft: Martijn van den Ende. Writing - Review & Editing: All authors. Supervision: Diane Rivet. Funding acquisition: Diane Rivet.

**Abstract** Owing to its deployment and sensing characteristics, Distributed Acoustic Sensing (DAS) has been touted as a promising technology for low-cost and low-latency Earthquake Early Warning (EEW). While preliminary experiments conducted by several research groups have yielded encouraging results, it must be acknowledged that these EEW feasibility studies were performed only on low-magnitude events. When exposed to the wavefield of a large magnitude earthquake (being the prime subject of EEW), the DAS strain rate recordings are likely to become highly distorted ("saturated") due to cycle skipping of the optical phase measurements, to an extent that the recorded data start to degrade to uniform random noise. This clearly poses a major challenge to EEW, as neither amplitude nor phase information can be readily extracted from saturated DAS data. In this study, we perform a detailed analysis of the dynamic range of DAS, both from theoretical and practical perspectives. We offer a set of criteria that need to be met for matching the DAS dynamic range with EEW targets, and we propose a computationally convenient method to quantify the information content of saturated recordings. We apply these methods to DAS data recorded offshore Chile, and identify several avenues for future research to improve the feasibility of DAS for EEW.

**Non-technical summary** Distributed Acoustic Sensing (DAS) is a relatively new technology that uses telecommunication (internet) cables to record vibrations in the ground. Because telecommunication cables are robust and available in many places, DAS could potentially be used for rapid detection of large earthquakes. Seismological institutes could use these data to send an alert for an imminent earthquake in near-real time, a concept that is known as Earthquake Early Warning (EEW). Unfortunately, there are some limitations of DAS that make it potentially unsuitable for EEW. The main limitation is that under strong ground shaking, the DAS measurements become unusable for the analysis of the earthquake. In this study we look into this limitation in detail, and we propose several mitigation strategies that could potentially make DAS more suitable for the analysis and EEW of large earthquakes.

**Resumen** La detección acústica distribuida (DAS) es una tecnología relativamente nueva que utiliza cables de telecomunicaciones (Internet) para registrar las vibraciones del suelo. Dado que los cables de telecomunicaciones son robustos y están disponibles en muchos lugares, los DAS podrían ser utilizados para detectar rápidamente grandes terremotos. Los institutos sismológicos podrían utilizar estos datos para enviar una alerta de terremoto inminente en tiempo casi real, un concepto que se conoce como Alerta Temprana de Terremoto (EEW). Desgraciadamente, el DAS tiene algunas limitaciones que lo hacen potencialmente inadecuado para la EEW. La principal es que, en caso de fuertes movimientos de suelo, las mediciones del DAS resultan inutilizables para el análisis del seísmo. En este estudio examinamos en detalle esta limitación y proponemos varias estrategias de mitigación que podrían hacer que el DAS fuera más adecuado para el análisis y la EEW de grandes terremotos.

## 1 Introduction

Earthquakes are among the most destructive of natural hazards, having claimed an estimated 2.5 million fatalities and over $900B in economic damages since 1900 (National Geophysical Data Center, 2023). Since short-term earthquake forecasting remains infeasible to this date, Earthquake Early Warning (EEW; Allen and Melgar, 2019) can be considered society's first line of defence against seismic hazard. The concept of EEW is based on the notion that telecommunication can trans-

mit information faster than the speed of seismic wave propagation. Seismic stations near the epicentre that experience strong ground shaking can trigger an alert that is transmitted to more distant localities, providing up to several seconds of lead time before the arrival of the seismic waves at those localities. During those precious few seconds, individuals can try to seek cover while various other mitigation measures can be (automatically) taken, such as the slowing down of trains and the initiation of emergency shutdown protocols of critical infrastructure. In order for EEW to be successful, a dense network of seismic sensors close to the epicentre

Production Editor:

Gareth Funning

Handling Editor:

Wenbo Wu

Bradley Lipovsky

Copy & Layout Editor:

Kirsty Bayliss

Signed reviewer(s):

Jimmy Atterholt

Qibin Shi

Received:

April 26, 2024

Accepted:

March 30, 2025

Published:

April 22, 2025

<footnote>*Corresponding author: martijn.vandenende@oca.eu</footnote>
is required. Unfortunately, achieving sufficient proximity can be challenging, especially for subduction zone settings that necessitate offshore sensor deployments to maximise the lead time.

In recent years, Distributed Acoustic Sensing (DAS; Hartog, 2017) has been considered as a potential solution for improving EEW coverage, particularly in subduction settings (Zhan, 2020; Farghal et al., 2022; Lior et al., 2023; Yin et al., 2023a). DAS is a fibre-optic sensing technique that converts a fibre-optic (e.g., telecom) cable into an array of equidistant vibration sensors by means of optical interferometry. One major advantage of DAS is that it can utilise existing telecom infrastructures (Lindsey et al., 2019; Sladen et al., 2019; Li et al., 2023a), drastically reducing the cost of deployment and maintenance of DAS-based seismic arrays, and granting access to environments that are inhospitable to conventional instrumentation (like offshore settings). However, even though some studies have already reported on the feasibility of DAS for earthquake seismology (Lior et al., 2023; Yin et al., 2023a), it is still a nascent technology of which the applicability and limitations need to be evaluated.

One limitation that is particularly relevant for EEW, is the response of the DAS instrument when subjected to strong ground motions. Conventional broadband seismometers exhibit a dynamic range and sensitivity that allow for the detection of low-amplitude ground motions down to the ambient noise level, but they saturate their measurements when the ground motions exceed a given range. Hence, EEW often relies on strong-motion accelerometric sensors that have lower sensitivity, but that do not saturate when subjected to strong motions. In principle, DAS exhibits a sensitivity that can be similar to that of conventional broadband seismometers (Lior et al., 2021), and it does not saturate its measurements in the same manner as seismometers. However, since DAS is an interferometric technique, it does suffer from data corruption due to failure to track the optical phase between consecutive laser pulses (also known as cycle skipping). Hence, while DAS does not saturate in the classical sense, it does exhibit data corruption when the optical fibre is exposed to sufficiently large strain rates (see Fig. 1).

While it is clear that the limited dynamic range of DAS poses a challenge for EEW applications, the issue itself has received relatively little attention. Diaz-Meza et al. (2023) reported on the observation of DAS “saturation” during tap tests, and proposed a detection/reconstruction algorithm that was able to recover from mild saturation in a specially engineered fibre (see Fig. 1d for an example that matches our definition of “mild”). Kong et al. (2022) proposed a Deep-Learning based solution for reliable phase unwrapping of single-channel Φ-OTDR DAS data, likewise applied to a (synthetic) scenario of mild saturation. Abukrat et al. (2023) remarked that saturation of near-source channels prevents accurate picking of seismic phases, and recommended the use of short fibres (i.e., short sensing distances) that permit the use of shorter gauge lengths (and correspondingly a larger dynamic range) while maintaining a reasonable data volume. A similar recommendation

was made by Viens et al. (2022), who were forced to exclude the strongly-saturated measurements of a nearby Mv 5.6 event from their analysis. However, reducing the sensing range of DAS would limit its usefulness in EEW, and should thus be avoided if possible.

The objective of the present study is to bring more attention to the limited dynamic range of DAS and its challenges for DAS-based EEW. We first conduct a theoretical analysis of the origin of DAS data saturation for a monochromatic oscillator, from which a criterion for the dynamic range emerges. We then extend this result to the case of a broadband earthquake spectrum, and test our model against the earthquake recordings generated by three 150-km long DAS cables that are located offshore Chile. Finally, we derive a metric for the degree of DAS data saturation that quantifies an upper saturation limit beyond which the data have become statistically indistinguishable from white noise. This metric is subsequently applied to the data of the 2023 Mww 6.6 Huasco earthquake. We conclude with recommendations for EEW-specific deployments and for future technological developments.

# 2 Analysis of DAS dynamic range

## 2.1 Simulating a DAS measurement

For a complete discussion of the dynamic range of DAS, we begin with a simple theoretical analysis. In what follows, it is helpful to recall the main operations by which a DAS interrogator converts a phase measurement into a measure of (local) strain rate (for a more in-depth discussion, see e.g. Hartog (2017); Lindsey et al. (2020)):

1. First, the interrogator sends a pulse into the sensing fibre and records the phase of the back-scattered light as a function of time. This so-called fast time axis can be converted into distance along the fibre using the speed of light in glass.

2. The pulsation is then repeated to obtain a subsequent phase measurement. The difference in phase at a given cable position is proportional to the length change of the fibre up to that position. In other words, by taking the time derivative along the time axis sampled by consecutive pulses (the “slow” time axis), a measure is obtained for the stretching rate of the fibre.

3. This stretching rate is a cumulative measurement up to a given point, and so to obtain a local measurement, the stretching rate is converted into strain rate by taking the spatial derivative (corresponding with the time derivative along the fast time axis).

For a quantitative analysis, the above procedure needs to be made more precise, which will be the objective for the remainder of this subsection.

Assuming an ideal instrument and perfectly linear response, a phase measurement Φ(t, x) corresponding to space-time lengthening d(t, x) of a fibre-optic cable is:

$$ \Phi(t, x) = \frac{4\pi\nu\xi}{\lambda}d(t, x) \quad (1) $$
<figure><img src="image_2.png" /><figcaption>**Figure 1** Experimental setting, and examples of DAS data saturation. (a) Overview of the ABYSS experiment, located along the Central Chilean margin. The route of each of the three DAS cables (CCN.N, SER.S, and SER.N) are indicated in black. The epicentre of the M<sub>ww</sub> 6.6 Huasco earthquake is marked by the orange star; (b) Time-series of the Huasco event recorded by a poorly-coupled DAS channel, which are not saturated; (c) Time-series of the Huasco event recorded by a well-coupled DAS channel. The data are strongly saturated, approaching the limit of uniform (white) noise; (d) Time-series of ocean gravity waves (swell) recorded close to shore. The data are mildly saturated, and could potentially be recovered by unwrapping; (e) Time-series of a vibrating segment, driven by ocean-bottom currents. Data saturation starts after 4 s. In panels (b)-(e), the dynamic range is indicated by the dotted lines.</figcaption></figure>

where $\nu$ denotes the refractive index of glass, $\xi$ the photo-elastic coefficient, and $\lambda$ the wavelength of light. Here, $t$ corresponds to the slow time axis. To facilitate arithmetic operations, from here on, the measurement is expressed in complex form: $m = e^{j\Phi}$ (with $j^2 \equiv -1$). In this formulation, it is guaranteed that $\Phi(t, x) = \arg\{m(t, x)\} \in [-\pi, +\pi)$ (with $\arg\{\cdot\}$ denoting the complex argument), eliminating the burden of explicitly considering the angle quadrant in the calculations. We will rely on this when we demonstrate a simple unwrapping algorithm (Appendix I).

The first processing step is to take the temporal derivative of the discrete phase measurement $\Phi(t_n, x_k)$:

$$ \begin{aligned} \dot{\Phi}(t_n, x_k) &\approx \frac{1}{\Delta t} [\Phi(t_n, x_k) - \Phi(t_{n-1}, x_k)] \quad (2) \\ &= \frac{1}{\Delta t} \arg \{m(t_n, x_k)m^*(t_{n-1}, x_k)\} \end{aligned} $$

with $t = [0, \Delta t, \dots, N\Delta t]^T$ and $x = [0, \Delta x, \dots, K\Delta x]^T$, and $m^*$ denoting the complex conjugate of $m$. As before, we define $\dot{m} = e^{j\Phi}$. To reduce the influence of measurement noise, $\dot{\Phi}(t_n, x_k)$ is typically averaged along the spatial dimension with a sliding window of fixed size, but

this is not an essential component of the analysis. Finally, the spatial gradient is taken to make the measurement local:

$$ \begin{aligned} \nabla \dot{\Phi}(t_n, x_k) &\approx \frac{1}{L\Delta x} [\dot{\Phi}(t_n, x_k) - \dot{\Phi}(t_n, x_{k-L})] \quad (3) \\ &= \frac{1}{L\Delta x} \arg \{\dot{m}(t_n, x_k)\dot{m}^*(t_n, x_{k-L})\} \end{aligned} $$

$\Delta x$ represents the spatial discretisation interval, which, depending on the specific interrogator model, may or may not correspond with the gauge length as defined by the optical pulse width. The integer $L > 0$ indicates the number of spatial samples over which the spatial gradient is computed, and effectively acts as an artificial gauge length. For the remainder of this work, we set $L = 1$. Correspondingly, the measurement of strain rate is obtained as:

$$ \begin{aligned} \dot{\varepsilon}(t_n, x_k) &= \frac{\lambda}{4\pi\nu\xi} \nabla \dot{\Phi}(t_n, x_k) \\ &= \frac{\lambda}{4\pi\nu\xi} \frac{\arg\{m_{n,k}m_{n-1,k}^*m_{n,k-1}^*m_{n-1,k-1}\}}{\Delta t \Delta x} \end{aligned} \quad (4) $$

with $m_{n,k} = m(t_n, x_k)$ for compactness of the notation.
## 2.2 The dynamic range of a monochromatic oscillator

Next, we derive expressions for the dynamic range of the measurement using a simple monochromatic oscillator as an example. Consider the following space-time dependence of the cable length $d(t, x)$:

$$d(t, x) = A \cos \left( 2\pi f \left[ t + \frac{x}{c} \right] \right) \quad (5)$$

with amplitude $A$, wave frequency $f$, and apparent phase velocity $c$. The exact expressions of the phase rate and its spatial gradient are:

$$\dot{\phi}(t, x) = -2\pi A f \frac{4\pi\nu\xi}{\lambda} \sin \left( 2\pi f \left[ t + \frac{x}{c} \right] \right) \quad (6a)$$

$$\nabla \dot{\phi}(t, x) = - \frac{4\pi^2 f^2 A}{c} \frac{4\pi\nu\xi}{\lambda} \cos \left( 2\pi f \left[ t + \frac{x}{c} \right] \right) \quad (6b)$$

The optical phase that defines the primary measurement of DAS is limited between $-\pi$ to $+\pi$, and consequently the difference between two consecutive phase measurements cannot exceed $\pi$. Hence, saturation of the DAS data can occur when $|\Phi(t_n, x_k) - \Phi(t_{n-1}, x_k)| > \pi$. In the case of a monochromatic oscillator, this occurs when:

$$\begin{aligned} |\Phi(t_n, x_k) - \Phi(t_{n-1}, x_k)| &\approx |\dot{\phi}(t_n, x_k)|\Delta t \\ &= 2\pi A f \frac{4\pi\nu\xi}{\lambda} \left|\sin\left(2\pi f \left[t + \frac{x}{c}\right]\right)\right| \Delta t > \pi \quad (7) \\ \Leftrightarrow \quad A &> \frac{1}{2\pi f} \frac{\lambda}{4\pi\nu\xi} \frac{\pi}{\Delta t} \frac{1}{\sin\left(2\pi f \left[t + \frac{x}{c}\right]\right)} \end{aligned}$$

Hence, the lower bound on the particle displacement amplitude that causes saturation of the time derivative is:

$$A_{\text{crit},t} = \frac{1}{2\pi f} \frac{\lambda}{4\pi\nu\xi} \frac{\pi}{\Delta t} \quad (8)$$

However, as can be seen in Fig. 2e, saturation of the temporal derivative does not necessarily lead to a saturation of the final measurement (strain rate). This is due to the spatial derivative, which does not yet experience similar saturation, and so even when the time derivative is saturated, its spatial derivative is not; in other words, the phase difference between $\dot{m}_{n,k}$ and $\dot{m}_{n,k-1}$ (in Eq. (3)) does not necessarily exceed $\pi$, even if the phase difference between $m_{n,k}$ and $m_{n-1,k}$ (in Eq. (2)) does. Continuing the same strategy adopted above, the saturation criterion for the monochromatic oscillator is given by:

$$\begin{aligned} |\dot{\Phi}(t_n, x_k) - \dot{\Phi}(t_n, x_{k-1})| \Delta t &\approx |\nabla \dot{\phi}(t_n, x_k)| \Delta x \Delta t \\ &= \frac{4\pi^2 f^2 A}{c} \frac{4\pi\nu\xi}{\lambda} \left|\cos\left(2\pi f \left[t + \frac{x}{c}\right]\right)\right| \Delta x \Delta t > \pi \quad (9) \\ \Leftrightarrow \quad A &> \frac{c}{4\pi^2 f^2} \frac{\lambda}{4\pi\nu\xi} \frac{\pi}{\Delta x \Delta t} \frac{1}{\left|\cos\left(2\pi f \left[t + \frac{x}{c}\right]\right)\right|} \end{aligned}$$

and:

$$A_{\text{crit},\hat{\epsilon}} = \underbrace{\frac{c}{4\pi^2 f^2}}_{\text{wavefield}} \frac{\lambda}{4\pi\nu\xi} \underbrace{\frac{\pi}{\Delta x \Delta t}}_{\text{instrument}} = \frac{c}{2\pi f \Delta x} A_{\text{crit},t} \quad (10)$$

In the above equations, the contributions of the wavefield (conversion from particle motion to strain), the

fibre (optical characteristics), and the instrument (interrogation settings) have been made explicit. Note that the first term labelled “wavefield” is specific to our choice of a monochromatic oscillator, but it can likewise be obtained through the well-known relationship between particle acceleration $a$ and strain rate $\dot{\epsilon}$, $a = c\dot{\epsilon}$ (Daley et al., 2016). In the spectral domain, the double integration that converts acceleration into displacement would correspond with multiplication of the strain rate spectrum with $(4\pi^2 f^2)^{-1}$, and so one would again obtain $c/4\pi^2 f^2$ as the conversion factor between peak ground displacement (A) and strain rate, where f would correspond with the characteristic frequency of a narrowband seismic source. However, the narrowband assumption is too restrictive for most seismic sources, and so the above criterion cannot directly be applied to earthquake data. Hence, we need to consider finite-source effects to obtain a criterion that is relevant in practice.

## 2.3 The dynamic range of a broadband signal

A conventional representation of an earthquake amplitude spectrum is given by the Brune spectrum, which, when accounting for frequency-dependent attenuation, reads (Brune, 1970; Anderson and Hough, 1984):

$$\Omega(f) = (2\pi f)^2 \frac{\Omega_0}{1 + \left[ \frac{f}{f_c} \right]^2} \exp(-\pi\alpha f) \quad (11)$$

where $\Omega$ denotes the acceleration spectrum (which is proportional to the strain rate spectrum) with reference spectrum $\Omega_0$, $f_c$ is the corner frequency, and $\alpha$ is an attenuation parameter. The reference spectrum and corner frequency are related to the seismic moment $M_0$ as (Madariaga, 1976; Shearer, 2011):

$$\Omega_0 = \frac{M_0 \Theta}{4\pi \rho c_s^3 R} \quad (12a)$$

$$f_c = k c_s \left( \frac{16}{7} \frac{\Delta \tau}{M_0} \right)^{1/3} \quad (12b)$$

In these expressions, $\rho$ denotes the mass density, $c_s$ the shear wave speed, $R$ the hypocentral distance, $k$ a geometric constant, and $\Delta\tau$ the mean stress drop across a circular crack. The parameter $\Theta$ comprises various contributions from the radiation pattern and free surface effects, including the broadside sensitivity (Martin et al., 2021) and cable-ground coupling in the case of DAS.

While the Brune acceleration spectrum gives an indication of which frequencies may exceed the saturation threshold, one must keep in mind that it is not the saturation of individual frequencies that is observed. In the time domain, it is the superposition of the contributions of each frequency that may ultimately exceed the dynamic range, and hence $\Omega(f)$ cannot be compared with $A_{\text{crit}}$ directly. A statistically robust alternative that translates the Brune spectrum into an equivalent time-domain signal amplitude is the root-mean squared acceleration $a_{\text{RMS}}$. Through the application of Parseval's theorem and subsequent simplification, Lior and Ziv
<figure><img src="image_3.png" /><figcaption>**Figure 2** Comparison between the optical phase induced by a monochromatic oscillator (orange line), and as measured by a theoretical DAS instrument (black line). The left column includes the phase (Φ), phase rate (Φ̇), and the gradient of the phase rate (∇Φ) for an oscillation amplitude of A = 1.5A crit,t < A crit,x. The right column includes the same quantities for A = 1.5A crit,x. Note that saturation of the final measurement (proportional to strain rate) only becomes saturated when A exceeds A crit,x, but not when it exceeds A crit,t.</figcaption></figure>

(2018) obtained the following approximation of a RMS for the attenuated Brune spectrum:

$$a_{\text{RMS}} = \frac{(2\pi f_c)^2 \Omega_0}{\sqrt{\pi\alpha T \left(1 + \left[\frac{2}{3}\right]^{\frac{1}{4}} \pi\alpha f_c\right)^2}} \quad (13)$$

Here, T denotes the duration of the time-domain signal for which a RMS is obtained. Using this expression, the (attenuated) spectral characteristics of the seismic source are translated into an equivalent signal amplitude a RMS that can be substituted into Eq. (10).

The next step is to compare a RMS with the phase saturation criterion. Note that Eq. (10) is defined in terms of particle displacement, whereas a RMS is a measure of acceleration. Hence, the “wavefield” term can be replaced with the apparent phase velocity c (since a = ċ); the same result is obtained by first differentiating Eq. (5) twice with respect to time and repeating the subsequent steps. This gives a crit = cλ (4νξΔxΔt)⁻¹. We then set the signal saturation threshold at a crit = 2a RMS, which is equivalent to having 95 % of the signal contained within the dynamic range (assuming a normal distribution of amplitudes), and express Ω₀ and f c in terms of M₀ (Eq. (12)). The expression that follows is cubic in M⁻¹, and so it permits an analytical solution of M₀ in terms of a RMS, R, etc., but the solution is too cumbersome to be of practical use. Instead, we recognise that Eq. (13) has two asymptotes around f c = [3/2]¹/₄ (πα)⁻¹, each of which permits a simple analytical solution of

M₀:

$$M_0 = \begin{cases} \left[ \frac{a_{\text{crit}} R \sqrt{\pi \alpha T}}{2(2\pi\sigma)^2\mu} \right]^3 & \text{for } f_c \ll \left[ \frac{3}{2} \right]^{\frac{1}{4}} (\pi\alpha)^{-1} \\ \alpha^2 \sqrt{\frac{2}{3}\pi\alpha T} \frac{R}{8\mu} a_{\text{crit}} & \text{for } f_c \gg \left[ \frac{3}{2} \right]^{\frac{1}{4}} (\pi\alpha)^{-1} \end{cases} \quad (14)$$

with:

$$\mu = \frac{\Theta}{4\pi\rho c_s^3} \\ \sigma = k c_s \left[ \frac{16}{7} \Delta\tau \right]^{\frac{1}{3}} \quad (15) \\ a_{\text{crit}} = \frac{c\lambda}{4\nu\xi\Delta x \Delta t}$$

The first asymptote represents the scaling of M₀ for large magnitude earthquakes (small f c), whereas the second represents that of small magnitude earthquakes. See Supplementary Figure S1 for a visualisation of these asymptotes as a function of M w, taking the parameters from Table 1.

Finally, following the conventional scaling between seismic moment M₀ and moment magnitude M w (Hanks and Kanamori, 1979), (14) can be used directly to compute the moment magnitude above which the DAS recordings become saturated:

$$M_w = \frac{2}{3} [\log_{10}(M_0) - 9.05] \quad (16)$$

Up to this point, the attenuation parameter α has been considered to be constant. However, conventionally
this parameter comprises the distance that a seismic ray has travelled through the attenuating medium in the form $\alpha = 2R(cQ)^{-1}$, where c and Q are the average phase speed and “quality factor” of the medium, respectively, and where the ray path is approximated by the hypocentral distance R.

To verify this relationship for the onset of DAS data saturation, we analyse the DAS data recorded for 73 events taken from the public catalogue of the *Centro Sismológico Nacional* (Universidad de Chile, 2012), ranging in magnitude from 2.5 up to 6.6. For each DAS channel we obtain the peak amplitude achieved within 30 s after the arrival of the first detectable phase arrival. To avoid strong ocean swell from affecting the results, we excluded the first 20 km of each cable. If the peak amplitude exceeds 90 % of the dynamic range, the given DAS channel is assumed to have been affected by saturation. Hence, this approach yields a binary classification of the onset of saturation as a function of catalogue magnitude and hypocentral distance. Other metrics, like those based on the total energy or the 90th-percentile, give very similar results as taking the maximum amplitude.

Aside from the parameters pertaining to the acquisition and fibre, we take the remaining parameters in Eq. (14) from Strumia et al. (2024) – see Table 1. A reasonable fit between the model and the observations (Fig. 3) is obtained when we assume an average coseismic stress drop of $\Delta\tau = 5$ MPa, which is a bit higher than what was inferred by Strumia et al. (2024). Then, there are numerous factors comprised in the “effective” radiation pattern $\Theta$, such as the orientation of the fault, the cable orientation, and the local velocity structure, all of which are unknown. Due to the broadside sensitivity of DAS, certain phases could potentially be recorded with almost zero amplitude, and therefore not trigger saturation of the data. However, since the analysis presented here considers 73 earthquakes distributed over a wide region, and three different cables with (someWHAT) variable geometry and distance to each seismic event, it is not physically realistic to assume a specific seismic phase, radiation pattern, or cable orientation. We therefore opt to consider a representative average of $\Theta$. Following the theoretical analysis of Strumia et al. (2024), the effective radiation pattern (including broadside sensitivity) averaged over the focal sphere and all possible cable orientations, takes a value of 0.2586 for the P-phase, and 0.2518 for the S-phase. Given the numerous simplifications and approximations made so far, we simply assume a value of $\Theta = 0.25$ to represent both phases.

A second justification for averaging the radiation pattern comes from the notion that the wavefield recorded by DAS is dominated by scattered phases. This is in part due to the higher sensitivity of DAS to slower phases, causing shallow scattered phases to be recorded with higher amplitude, and in part due to the shallow sedimentary structure that is typical for marine environments, promoting scattering of incoming seismic waves (Trabattoni et al., 2024). By taking a time window of a certain duration (e.g. 10 or 30 seconds), it is likely that the recorded wavefield will comprise many scattered

<table><caption>**Table 1** Selected parameters for Eq. (17) and Fig. 3</caption><thead><tr><th>Quantity</th><th>Symbol</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>App. wave speed</td><td>c</td><td>400</td><td>m s<sup>-1</sup></td></tr><tr><td>S-wave speed</td><td>c<sub>s</sub></td><td>2500</td><td>m s<sup>-1</sup></td></tr><tr><td>Geometric factor</td><td>k</td><td>0.26</td><td>-</td></tr><tr><td>Stress drop</td><td>Δτ</td><td>5 × 10<sup>6</sup></td><td>Pa</td></tr><tr><td>Radiation pattern</td><td>Θ</td><td>0.25</td><td>-</td></tr><tr><td>Quality factor</td><td>Q</td><td>800</td><td>-</td></tr><tr><td>Mass density</td><td>ρ</td><td>2700</td><td>km<sup>-3</sup></td></tr><tr><td>Optical wavelength</td><td>λ</td><td>1550 × 10<sup>-9</sup></td><td>m</td></tr><tr><td>Refractive index</td><td>ν</td><td>1.44</td><td>-</td></tr><tr><td>Photo-elastic coeff.</td><td>ξ</td><td>0.79</td><td>rad</td></tr><tr><td>Gauge length</td><td>Δx</td><td>30</td><td>m</td></tr><tr><td>Time sampling rate</td><td>Δt</td><td>16 × 10<sup>-3</sup></td><td>s</td></tr><tr><td>Time window</td><td>T</td><td>30</td><td>s</td></tr></tbody></table>

arrivals with a relatively slow apparent velocity. If the scatterers are assumed to be isotropically distributed, the effective $\Theta$ will again take a value that is an average over many different orientations. Moreover, this would constrain the apparent wave speed c to be representative for the shallow sedimentary structures underlying the cable; here, we take a value of 400 m s$^{-1}$.

Considering that not all events conform to representative average parameters that are assumed for the model, deviations from the predicted saturation threshold are expected for individual events. What is particularly clear from Fig. 3 is that the scaling of this threshold is sensitive to the attenuation: by setting a fixed $\alpha = 0.1$ s, a distance scaling is obtained that does not match the observations, underestimating the saturation potential of proximal events and overestimating that of distant events. By accounting for the length of the ray path, a much more reasonable scaling is obtained. For the attenuation parameters that describe the data well with $\alpha = f(R)$, it is found that $f_c \ll [\frac{3}{2}]^{\frac{1}{4}} (\pi\alpha)^{-1}$ for all events, and so only the first asymptote of Eq. (14) is practically relevant. The analysis of this section can thus be summarised with the following expression:

$$M_w = 2 \log_{10} \underbrace{\left(a_{\text{crit}} \sqrt{\frac{T}{2\pi}}\right)}_{\text{acquisition}} - 4 \log_{10} \underbrace{\left(k \left[\frac{16}{7} \Delta\tau\right]^{\frac{1}{3}} \sqrt{\Theta}\right)}_{\text{source}} + \log_{10} \underbrace{\left(\frac{\rho^2 c_s}{Q}\right)}_{\text{medium}} + 3 \log_{10}(R) - \frac{2}{3} \cdot 9.05 \quad (17)$$

The scaling of this final expression is consistent with the data in Fig. 3b, displaying a constant magnitude-distance scaling consistent with this result ($M_0 \propto R^3$).

One final observation is that for channels that are in the proximity (< 50 km) of the hypocentre, saturation is observed (and predicted) to occur for magnitudes as low as 3. This clearly presents a challenge to EEW efforts, which are most effective when deployed in the vicinity of the seismic source. However, Fig. 3 only considers the onset of saturation, at which point the recordings retain much of their original information content (i.e., the case of “mild” saturation as shown in Fig. 1d). In
<figure><img src="image_4.png" /><figcaption>**Figure 3** Observed and predicted DAS data saturation induced by 73 seismic events. Each DAS channel is classified as “saturated” (orange) or “not saturated” (green) based on the recorded peak amplitude relative to its dynamic range. Eq. (14) is plotted taking a fixed $\alpha = 0.1$ s (dotted line) or as a function of hypocentral distance $\alpha = 2R(cQ)^{-1}$ (dashed line). The transition to white noise as described in Section 3 is given by the solid black line. Panel (b) displays the same information as panel (a), but on a logarithmic distance scale to highlight the more proximal, smaller magnitude events.</figcaption></figure>

the next section, we will consider how rapidly this information content degrades with increasing amplitude. Moreover, we define a metric for the degree of saturation, and we use this metric to estimate an upper bound beyond which no useful information can be (theoretically) recovered.

## 3 Spectral distortion and saturation metric

As discussed in the previous section, signals that exceed the DAS dynamic range cause distortions in the recordings, which negatively affects their useful information content. However, this transformation is not instantaneous: while strongly saturated data (Fig. 1c) bear little resemblance with the true underlying signal, weakly saturated data (Fig. 1b) could potentially be unwrapped or have their spectra analysed to extract the corner frequency (Strumia et al., 2024). To see the effect of phase saturation on a broadband earthquake spectrum, we use Eq. (11) to generate a synthetic source spectrum with uniform random phase, and convert it into a time series with an inverse Fourier transformation. We then scale the signal peak amplitude by a factor in the range of 1.5 to 5.0 times the dynamic range, and wrap the signal. The spectra that are observed after this synthetic saturation are shown in Fig. 4. For mild saturation (Fig. 4a), the spectral information around the corner frequency remains mostly unaffected, but a prominent low-frequency plateau emerges. As the synthetic time series becomes increasingly more saturated (panels b and c), this plateau increases in amplitude while the amplitude of the spectral peak diminishes, until the entire spectrum starts to approach a white noise spectrum (panels d and e).

To get a better grasp on this spectral behaviour, we first note that any saturated time series $\overline{m}(t)$ with a dynamic range of $\pm A$ can be decomposed into a superposition of the original signal $m(t)$, and a set of rectangle functions of amplitude $\pm 2A$ and width $w$, i.e.:

$$ \overline{m}(t) = m(t) + 2A \sum_{i=0}^{I} \mathrm{rect}_{w_i}(t - t_i) - 2A \sum_{j=0}^{J} \mathrm{rect}_{w_j}(t - t_j) \quad (18) $$

with

$$ \mathrm{rect}_w(t - t_0) = H \left\{ t - \left(t_0 - \frac{w}{2}\right) \right\} - H \left\{ t - \left(t_0 + \frac{w}{2}\right) \right\} \quad (19) $$

and $H\{\cdot\}$ denoting the Heaviside step-function. As a consequence, the spectrum of $\overline{m}(t)$ can be expressed as the weighted summation of $\Omega(f)$ (i.e., the spectrum of the true signal) and the spectrum corresponding to the superposition of rectangle functions $\Lambda(f)$:

$$ \begin{aligned} |\mathcal{F}(\overline{m})|^2 &= (\Omega + 2A\Lambda)(\Omega + 2A\Lambda)^* \\ &= |\Omega|^2 + 4A^2|\Lambda|^2 + 4A\mathrm{Re}(\Omega\Lambda^*) \end{aligned} \quad (20) $$

Here, $\mathcal{F}(\cdot)$ is the Fourier transform. The last term on the right-hand side denotes the real component of the cross-correlation between the true signal and the summation of rectangle functions, which in the spectral domain can be expressed as the multiplication of $\Omega$ with the complex conjugate of $\Lambda$. For mild saturation, cycle skipping is rare and only a few rectangle functions are needed to satisfy Eq. (18). The spectrum of a single rectangle function is given by the cardinal sine (or sinc) function, i.e. $|\mathcal{F}(\mathrm{rect}_w)| = |\sin(\pi fw^{-1}) / \pi f|$, and so $\Lambda(f)$ is well represented by the envelope of a sinc function. We show this in Fig. 4a, where the sinc spectrum is overlain on the spectrum of rectangle functions.
<figure><img src="image_5.png" /><figcaption>**Figure 4** (a)-(e) Synthetic acceleration spectra that undergo progressively more distortion, with peak amplitudes that range from 1.5 to 5.0 times the dynamic range. (f)-(j) The cumulative distribution function (CDF) of the synthetic signal amplitudes (normalised by the dynamic range), corresponding with the spectra shown in (a)-(e). As the signals become increasingly saturated, they gradually approach a uniform distribution (indicated by the dotted line). The saturation metric $D_T$ is as denoted in each panel.</figcaption></figure>

As given by Eq. (18), the low-frequency plateau, as well as the high-frequency distortion of the observed amplitude spectrum, originates from the addition of the sinc spectrum.

As the degree of saturation increases, the superposition of rectangle functions is no longer quasi-random, as it needs to cancel out the true earthquake spectrum to produce the ultimately observed, saturated spectrum. In other words, the contribution of Re($\Omega\Lambda^*$) becomes significant and $\Lambda$ is no longer described by a sinc function (Fig. 4b-e), resulting in a non-trivial superposition of spectra. However, even though the observed signal spectrum is no longer recognisable as a Brune spectrum, it has not yet become completely uniform (Fig. 4d-e). To quantify this, we consider how the amplitude distribution of $\bar{m}(t)$ approaches that of a uniform distribution (see Fig. 4f-j). The cumulative distribution function (CDF) of a uniform distribution $U[0, 1]$ is simply given by $C_U(x) = x (0 \le x < 1)$, and so when $\bar{m}$ is suitably scaled between [0, 1], the distance between the observed CDF ($C(x)$) and the uniform CDF can be conveniently defined as:

$$D = -\log_{10} \left( 3 \int_{0}^{1} [C(x) - x]^2 dx \right) \quad (21)$$

When $C(x)$ is perfectly non-uniform, the integral evaluates

to $1/3$, so that $D = 0$. As $\bar{m}$ approaches a uniform distribution, the integral evaluates to zero, and $D \to \infty$. In practice, the finite precision of the empirical CDF limits how close $C(x)$ can approach $C_U$, which in itself is a function of the number of samples $N_T$ contained within a given time window (in other words, the integration spacing $dx$). We find that the expectation of the upper limit of $D$ scales as $\log_{10}(3N_T)$, and so we define $D_T = D \log_{10}(3N_T)$, such that $0 < E[D_T] \le 1$. This renders the saturation metric independent of the arbitrary choice of $N_T$. The procedure of quantifying the degree of saturation of a given DAS channel is then to:

1. select a fixed time window (e.g., 30 s after the P-arrival),

2. compute the empirical CDF of the time series,

3. scale the data by the dynamic range such that all values fall between 0 and 1,

4. evaluate Eq. (21) and scale by $\log_{10}(3N_T)$.

Now that a suitable metric for the degree of saturation has been defined, we evaluate to what extent this metric can be applied in practice. We select an earthquake recorded by the SER.N cable with good signal-to-noise ratio, but without causing the data to saturate. For a selected number of DAS channels with good coupling and
<figure><img src="image_6.png" /><figcaption>**Figure 5** (a) Progressive degradation of an earthquake waveform towards white noise, with corresponding *D*<sub>T</sub>-values. (b) The evolution of the saturation metric (*D*<sub>T</sub>) with increasing signal peak amplitude (before wrapping). *D*<sub>T</sub> plateaus at an amplitude factor of 20, meaning that the saturated waveforms become statistically indistinguishable from white noise.</figcaption></figure>

low ambient noise levels, we normalise the data by the 90th-percentile value of each DAS channel. We then artificially saturate the data by scaling with a factor A and wrapping around the dynamic range, and measure D<sub>T</sub> as a function of A. Due to variations in the signal characteristics of each channel, the measured D<sub>T</sub>-values vary from one channel to the next. This allows us to estimate the expected variation of D<sub>T</sub> for a given earthquake, and consequently the precision of estimating A from D<sub>T</sub>. As seen in Fig. 5, the transition towards white noise is completed at around an amplitude factor of 20, beyond which D<sub>T</sub> plateaus. The exact point at which this occurs depends on the data window selection, as the S-wave train saturates more quickly than the (lower amplitude) P-wave train. Taking a factor 20 as a reasonable estimate, it is expected that the DAS data become “fully saturated” within 0.9 magnitude units (=$\frac{2}{3} \log_{10} 20$) above the initial saturation threshold. For reference, this upper bound on the saturation is included in Fig. 3.

To verify the existence of the upper bound in our data, we evaluate *D*<sub>*T*</sub> for the M<sub>ww</sub> 6.6 Huasco event, which is the largest magnitude event in the data set. The SER.N and SER.S cables were located within a hypocentral distance of less than 200 km, while the CCN.N cable was located around a distance of 400 km. By comparison with the fully-saturated bound in Fig. 3, one can see that the distal CCN.N cable is close to the predicted transition to becoming fully saturated, whereas the SER.N and SER.S cables are expected to have become fully saturated. We estimate *D*<sub>*T*</sub> for each channel individually, and observe how *D*<sub>*T*</sub> changes as a function of hypocentral distance and time (Fig. 6). When estimating *D*<sub>*T*</sub> over a time window spanning 10 s after the P-wave arrival at each channel, we observe that the CCN.N cable at around 400 km hypocentral distance is tightly clustered around *D*<sub>*T*</sub> ≈ 0.25, indicating that no (or very little) data saturation occurs. By contrast, the SER.N and SER.S cables (up to 200 km distance) exhibit *D*<sub>*T*</sub>-values that exceed 0.6, indicating that the data are saturated, though potentially retain some useful information. However, when extending the time window to 30 s and 60 s after the P-arrival, the SER.N and SER.S data approach the

white noise limit at *D*<sub>*T*</sub> ≈ 1. Taking the same time window, the CCN.N data start to exhibit a *D*<sub>*T*</sub>-value greater than 0.25, a trend that continues when extending the time window up to 60 s, reaching up to *D*<sub>*T*</sub> = 0.6.

These observations confirm that the plateau of *D*<sub>*T*</sub> observed in the quasi-synthetic analysis (Fig. 5b) is a phenomenon that manifests itself when the cable is subjected to strong ground motions, resulting in DAS recordings that have been fully reduced to white noise. The magnitude-distance criterion that marks the fully-saturated transition (Fig. 3) seems somewhat too stringent since the CCN.N cable, which was predicted to become fully saturated, exhibits intermediate *D*<sub>*T*</sub>-values up to 0.6. Based on Fig. 5b, in order to reach *D*<sub>*T*</sub> ≥ 0.9 the ground acceleration at CCN.N would have needed to be stronger by about a factor 3 (equivalent to a magnitude increase of 0.3). However, given the numerous simplifications made up to this point, and the uncertainties in the parameters that enter Eq. (14), we believe that the predictive power of the proposed saturation criteria is acceptable.

With regard to Fig. 6, we make one final observation, being that the lowest measured *D*<sub>*T*</sub>-values remain around 0.25 even for the SER.N and SER.S cables that reach the white noise transition. These are sections of the DAS cable that are poorly coupled, and that record little to no strain induced by the body waves. Some of these channels record only a fraction of the seismic energy, such that they remain unsaturated or become only mildly saturated, if at all (see Fig. 1b). These channels could be used to recover information that was lost by fully-saturated channels, and possibly play a critical role in EEW and near-field analysis of large seismic events. In this respect, well-coupled and poorly-coupled DAS segments would play a similar role as broadband and strong-motion sensors in conventional seismic networks. However, whether the data recorded by poorly coupled sections provides a faithful representation of the seismic wavefield (up to a scaling factor) still needs to be investigated.
<figure><img src="image_7.png" /><figcaption>**Figure 6** Saturation metric for the M<sub>ww</sub> 6.6 Huasco event. In each panel, *D*<sub>T</sub> is computed as a function of hypocentral distance, over a time window that extends 10, 30, or 60 seconds from the first arrival at each position. The colour intensity of the hexagonal bins is proportional to the number of data points in each bin. The onset of saturation at *D*<sub>T</sub> = 0.25 and the plateau at *D*<sub>T</sub> = 1 are indicated by the dotted and dashed lines, respectively.</figcaption></figure>

# 4 Implications for EEW

From the theoretical analysis laid out in Section 2, as well as from the empirical observations made in the previous section, it becomes clear that the limited dynamic range of DAS data presents a major obstacle for the use of DAS for EEW. Even though DAS has the potential to provide near-source instrumentation, which would maximise the warning time given by EEW (Lior et al., 2023), rapid saturation of the recordings prevents the accurate extraction of amplitude information that underlies local magnitude scaling relationships (Yin et al., 2023b). We found that for hypocentral distances of around 20 km, earthquakes of a magnitude as low as 2.5 could cause data saturation, which underlines the existing challenges that DAS-based EEW still needs to overcome.

Fortunately, there are several perspectives that could redeem DAS as an effective method for EEW:

* Firstly, while the amplitude information of the DAS recordings may not be informative, the P-wave onset times can be used to obtain a rapid first estimation of the hypocentre (Yin et al., 2023a). Given a known seismic source location, it is possible to estimate the source magnitude from the recordings of a single strong-motion sensor. This alleviates the need for several strong-motion sensors to be triggered, which is a requirement to obtain a source location (and corresponding magnitude) in conventional EEW systems. As a result, an alert could be issued as soon as the nearest strong-motion sensor exceeds a trigger limit, reducing the system latency by several seconds (Minson et al., 2018; Peng et al., 2021).

* Owing to imperfect coupling between the cable and the ground, DAS data commonly exhibit large variations in the spatial distribution of seismic wave amplitudes. While these poorly-coupled sections are typically considered a nuisance for weak-motion analyses, they may be vital to estimating ground

motion amplitudes that may otherwise have caused the data to become fully saturated. However, it needs to be investigated whether the response of these sections is linear (i.e., there exists a constant scaling between the imposed and recorded amplitudes), and whether they interact with soil nonlinearities (Viens et al., 2022).

* Particularly for submarine DAS, the shallow, unconsolidated sediment cover can cause major distortion of the seismic wavefield, including phase splitting (Trabattoni et al., 2024). As observed in a previous DAS study offshore Chile (Trabattoni et al., 2023), the arrival of the faint direct P-wave is commonly obscured by subsequent converted arrivals. This is a consequence of the measurement principle of DAS, measuring strains instead of particle motions, but it can be corrected by performing a spatial integration. After spatial integration, the direct P-wave can be distinguished from the P-to-S conversion that usually follows within 1 s. While the very first second of an earthquake trace is insufficient to establish a magnitude (Meier et al., 2017), the polarity and amplitude of the first motion can be used to construct an initial focal mechanism (Li et al., 2023b). In turn, this focal mechanism may serve to distinguish megathrust events occurring on the plate interface from those occurring in the overriding accretionary wedge or in the outer rise. It has been observed that such intraplate events can trigger disproportionally large tsunamis (Cummins and Kaneda, 2000; Hananto et al., 2020), and so DAS may contribute to improve tsunami early warnings.

Beside these future avenues for exploration, we offer the following recommendations for DAS-based EEW in its current form:

* As is apparent from Eq. (10), and as was recommended by previous studies (Viens et al., 2022; Abukrat et al., 2023), increasing the spatio-temporal sampling rate of the acquisition increases
the saturation amplitude proportionally. As it may not be feasible to record and store the acquired data without downsampling (owing to limited storage capacity or bandwidth), operating in trigger-mode could be a viable solution. However, one must keep in mind that the logarithmic magnitude scaling implies severely diminishing returns: by simultaneously decreasing the gauge length and temporal sampling period by a factor 10 (i.e, increasing the data density by a factor 100), the magnitude threshold is raised by only 1.3 magnitude units. Moreover, there are practical limitations to the data resolution which stem from the laser pulse duration and the total length of the sensed fibre. These prevent one from setting an arbitrary acquisition sampling rate in an attempt to extend the saturation threshold to the desired earthquake magnitude range. Another point of consideration is that for a fixed bit-resolution of the data (typically stored as 16 or 32-bit integers), the sensitivity scales proportionally with the dynamic range. For DAS arrays that serve both EEW and general seismological purposes (microseismicity monitoring, ambient noise correlation, etc.), a fixed bit-resolution may render the data unsuitable for either application. As an example, consider a dynamic range that permits unsaturated recordings of an $M_w$ 7 event at 50 km distance. A 16-bit discretisation (values in the range $\pm 2^{15}$) then implies that events of $M_w$ 4 or lower will generate strains that fall within the bit-precision, and are therefore rendered undetectable.

* As we have shown in Section 3, it is possible to make an estimation of the peak ground motion amplitude even when it exceeds the dynamic range. As long as the data distribution is significantly different from that of a uniform distribution, one could map the observed saturation metric ($D_T$) back into a peak amplitude following the relationship depicted in Fig. 5b. By doing so, the dynamic range can be artificially extended by an amount that corresponds with almost one unit of $M_w$. If the objective is to use an amplitude-based magnitude-distance relationship, extending the dynamic range also implies that reliable estimates of magnitude can be made much closer to the seismic source, which translates into faster alert times.

* By implementing (real-time) unwrapping algorithms, recordings that suffer from saturation could be fully restored. For mildly saturated data, one could attempt to detect sudden jumps in the data and add/subtract the dynamic range value to unwrap the saturated recordings (e.g. Diaz-Meza et al., 2023). However, this algorithm is severely limited and prone to errors once the true ground motion amplitudes far exceed the dynamic range. Instead, one should opt for a gradient-based unwrapping method, which tracks and integrates the phase-gradient of the recordings. We present a simple example of such an algorithm in Appendix I. More sophisticated gradient-tracking algorithms could greatly improve upon the unwrapping performance,

potentially increasing the magnitude detection threshold by more than one magnitude unit. The unwrapping should be performed prior to downsampling, as to retain the maximum sampling density for the gradient estimation, and not be affected by the step-response of anti-aliasing filters.

* One of the appealing features of DAS, is that it can leverage existing telecom fibre networks. However, for the specific purpose of EEW, it may be worthwhile (or necessary) to design and deploy customised cables that exhibit a lower sensitivity. This can be achieved either through changes in the optical characteristics (increasing optical wavelength or decreasing refractive index), or through geometry (e.g., helically-wound cables with an optimised pitch that can accommodate axial deformation with a smaller change in optical path length).

## 5 Conclusions

In this study, we presented an in-depth analysis of the dynamic range of Distributed Acoustic Sensing (DAS), with a primary focus on Earthquake Early Warning (EEW). While DAS offers many advantages over conventional instrumentation used by EEW systems, it suffers from data degradation (“saturation”) when the ground motions exceed its dynamic range. We derived several criteria that describe the tendency of DAS data saturation as a function of seismic moment magnitude and hypocentral distance, verified with empirical observations. From this we conclude that for typical DAS acquisition settings (gauge length and sampling rate), saturation occurs for a range of moment magnitudes and distances that is overly restrictive, indicating that present-day DAS technologies may not be suitable for EEW purposes. Furthermore, we proposed a metric for the degree of saturation, that may help to artificially extend the dynamic range of DAS by a factor equivalent to 0.9 units of $M_w$. Even though this is a significant improvement, it remains insufficient for the near-field analysis of earthquake magnitudes that are of interest of EEW (typically $M_w$ 6 or higher).

These observations indicate that technological advances still need to be made before DAS could replace conventional strong-motion instrumentation in EEW systems. Nonetheless, DAS can provide complementary information that helps to rapidly establish the seismic source location and focal mechanism, and dedicated signal processing techniques and custom cable designs could remedy the limited dynamic range, allowing DAS to contribute faithful amplitude and phase information to benefit EEW.

## Acknowledgements

We thank the editor, Bradley Lipovsky, and two anonymous reviewers for their suggestions to help improve the manuscript. This work was supported by the European Research Council (ERC) under the European
Union's Horizon 2020 research and innovation programme (grant agreement No. 101041092 – ABYSS). We acknowledge the assistance of Alcatel Submarine Networks (ASN) and GTD Telecom, as well as the Centro Sismológico Nacional (CSN). The following Python libraries were used for processing and visualisation of the data: JAX (Bradbury et al., 2020), NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2019), Pandas (Pandas Development Team, 2020), Cartopy (Office, 2015), Matplotlib (Hunter, 2007), Scientific Colourmaps (Crameri et al., 2020).

## Data and code availability

The data and scripts needed to reproduce the figures in this manuscript can be found in the following repository: https://doi.org/10.5281/zenodo.10993306

## Competing interests

The authors declare no competing interests.

## References

Abukrat, Y., Sinitsyn, P., Reshef, M., and Lellouch, A. Applications and Limitations of Distributed Acoustic Sensing in Shallow Seismic Surveys and Monitoring. *GEOPHYSICS*, 88(6):WC1-WC12, Nov. 2023. doi: 10.1190/geo2022-0574.1.

Allen, R. M. and Melgar, D. Earthquake Early Warning: Advances, Scientific Challenges, and Societal Needs. *Annual Review of Earth and Planetary Sciences*, 47(1):361–388, 2019. doi: 10.1146/annurev-earth-053018-060457.

Anderson, J. G. and Hough, S. E. A Model for the Shape of the Fourier Amplitude Spectrum of Acceleration at High Frequencies. *Bulletin of the Seismological Society of America*, 74(5): 1969–1993, Oct. 1984. doi: 10.1785/BSSA0740051969.

Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., and Zhang, Q. JAX: Composable Transformations of Python+NumPy Programs, 2020. http://github.com/jax-ml/jax.

Brune, J. N. Tectonic Stress and the Spectra of Seismic Shear Waves from Earthquakes. *Journal of Geophysical Research (1896-1977)*, 75(26):4997–5009, 1970. doi: 10.1029/JB075i026p04997.

Crameri, F., Shephard, G. E., and Heron, P. J. The Misuse of Colour in Science Communication. *Nature Communications*, 11(1): 5444, Oct. 2020. doi: 10.1038/s41467-020-19160-7.

Cummins, P. R. and Kaneda, Y. Possible Splay Fault Slip during the 1946 Nankai Earthquake. *Geophysical Research Letters*, 27(17): 2725–2728, 2000. doi: 10.1029/1999GL011139.

Daley, T. M., Miller, D. E., Dodds, K., Cook, P., and Freifeld, B. M. Field Testing of Modular Borehole Monitoring with Simultaneous Distributed Acoustic Sensing and Geophone Vertical Seismic Profiles at Citronelle, Alabama. *Geophysical Prospecting*, 64 (5):1318–1334, 2016. doi: 10.1111/1365-2478.12324.

Diaz-Meza, S., Jousset, P., Currenti, G., Wollin, C., Krawczyk, C., Clarke, A., and Chalari, A. On the Comparison of Records from Standard and Engineered Fiber Optic Cables at Etna Volcano (Italy). *Sensors*, 23(7):3735, Jan. 2023. doi: 10.3390/s23073735.

Farghal, N. S., Saunders, J. K., and Parker, G. A. The Potential of Using Fiber Optic Distributed Acoustic Sensing (DAS) in

Earthquake Early Warning Applications. *Bulletin of the Seismological Society of America*, 112(3):1416-1435, Apr. 2022. doi: 10.1785/0120210214.

Hananto, N. D., Leclerc, F., Li, L., Etchebes, M., Carton, H., Tapponnier, P., Qin, Y., Avianto, P., Singh, S. C., and Wei, S. Tsunami Earthquakes: Vertical Pop-up Expulsion at the Forefront of Subduction Megathrust. *Earth and Planetary Science Letters*, 538: 116197, May 2020. doi: 10.1016/j.epsl.2020.116197.

Hanks, T. C. and Kanamori, H. A Moment Magnitude Scale. *Journal of Geophysical Research: Solid Earth*, 84(B5):2348-2350, 1979. doi: 10.1029/JB084iB05p02348.

Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J. F., Wiebe, M., Peterson, P., Gérard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., and Oliphant, T. E. Array Programming with NumPy. *Nature*, 585(7825):357-362, Sept. 2020. doi: 10.1038/s41586-020-2649-2.

Hartog, A. H. *An Introduction to Distributed Optical Fibre Sensors*. CRC Press, May 2017. doi: 10.1201/9781315119014.

Hunter, J. D. Matplotlib: A 2D Graphics Environment. *Computing in Science & Engineering*, 9(3):90-95, 2007. doi: 10.1109/M-CSE.2007.55.

Kong, L., Pan, X., Ren, Z., and Cui, K. Robust One-Dimensional Phase Unwrapping Algorithm Based on LSTM Network With Reduced Parameter Number. *Journal of Lightwave Technology*, 40 (19):6560-6567, Oct. 2022. doi: 10.1109/JLT.2022.3195932.

Li, J., Kim, T., Lapusta, N., Biondi, E., and Zhan, Z. The Break of Earthquake Asperities Imaged by Distributed Acoustic Sensing. *Nature*, 620(7975):800-806, Aug. 2023a. doi: 10.1038/s41586-023-06227-w.

Li, J., Zhu, W., Biondi, E., and Zhan, Z. Earthquake Focal Mechanisms with Distributed Acoustic Sensing. *Nature Communications*, 14(1):4181, July 2023b. doi: 10.1038/s41467-023-39639-3.

Lindsey, N. J., Dawe, T. C., and Ajo-Franklin, J. B. Illuminating Seafloor Faults and Ocean Dynamics with Dark Fiber Distributed Acoustic Sensing. *Science*, 366(6469):1103-1107, Nov. 2019. doi: 10.1126/science.aay5881.

Lindsey, N. J., Rademacher, H., and Ajo-Franklin, J. B. On the Broadband Instrument Response of Fiber-Optic DAS Arrays. *Journal of Geophysical Research: Solid Earth*, 125(2): e2019JB018145, 2020. doi: 10.1029/2019JB018145.

Lior, I. and Ziv, A. The Relation Between Ground Motion, Earthquake Source Parameters, and Attenuation: Implications for Source Parameter Inversion and Ground Motion Prediction Equations. *Journal of Geophysical Research: Solid Earth*, 123(7): 5886-5901, 2018. doi: 10.1029/2018JB015504.

Lior, I., Sladen, A., Rivet, D., Ampuero, J.-P., Hello, Y., Becerril, C., Martins, H. F., Lamare, P., Jestin, C., Tsagkli, S., and Markou, C. On the Detection Capabilities of Underwater DAS. *Journal of Geophysical Research: Solid Earth*, n/a(n/a):e2020JB020925, 2021. doi: 10.1029/2020JB020925.

Lior, I., Rivet, D., Ampuero, J.-P., Sladen, A., Barrientos, S., Sánchez-Olavarría, R., Villarroel Opazo, G. A., and Bustamante Prado, J. A. Magnitude Estimation and Ground Motion Prediction to Harness Fiber Optic Distributed Acoustic Sensing for Earthquake Early Warning. *Scientific Reports*, 13(1):424, Jan. 2023. doi: 10.1038/s41598-023-27444-3.

Madariaga, R. Dynamics of an Expanding Circular Fault. *Bull. Seismol. Soc. Am*, pages 639-666, 1976.

Martin, E. R., Lindsey, N. J., Ajo-Franklin, J. B., and Biondi, B. L. Introduction to Interferometry of Fiber-Optic Strain Measurement.
ments. In Li, Y., Karrenbach, M., and Ajo-Franklin, J. B., editors, *Geophysical Monograph Series*, pages 111–129. Wiley, 1 edition, Dec. 2021. doi: 10.1002/9781119521808.ch9.

Meier, M.-A., Ampuero, J. P., and Heaton, T. H. The Hidden Simplicity of Subduction Megathrust Earthquakes. *Science*, 357(6357): 1277–1281, Sept. 2017. doi: 10.1126/science.aan5643.

Minson, S. E., Meier, M.-A., Baltay, A. S., Hanks, T. C., and Cochran, E. S. The Limits of Earthquake Early Warning: Timeliness of Ground Motion Estimates. *Science Advances*, 4(3):eaaq0504, Mar. 2018. doi: 10.1126/sciadv.aaq0504.

National Geophysical Data Center. Global Significant Earthquake Database, 2023. doi: 10.7289/V5TD9V7K.

Office, M. Cartopy: A Cartographic Python Library with a Matplotlib Interface, 2015. https://scitools.org.uk/cartopy.

Pandas Development Team, t. Pandas-Dev/Pandas: Pandas. Zenodo, 2020. doi: 10.5281/ZENODO.3509134.

Peng, C., Jiang, P., Ma, Q., Wu, P., Su, J., Zheng, Y., and Yang, J. Performance Evaluation of an Earthquake Early Warning System in the 2019–2020 M6.0 Changning, Sichuan, China, Seismic Sequence. *Frontiers in Earth Science*, 9, July 2021. doi: 10.3389/feart.2021.699941.

Shearer, P. M. *Introduction to Seismology*. Cambridge Univ. Press, Cambridge, 2. ed., repr. with corr edition, 2011.

Sladen, A., Rivet, D., Ampuero, J. P., De Barros, L., Hello, Y., Calbris, G., and Lamare, P. Distributed Sensing of Earthquakes and Ocean-Solid Earth Interactions on Seafloor Telecom Cables. *Nature Communications*, 10(1):1–8, Dec. 2019. doi: 10.1038/s41467-019-13793-z.

Strumia, C., Trabattoni, A., Supino, M., Baillet, M., Rivet, D., and Festa, G. Sensing Optical Fibers for Earthquake Source Characterization Using Raw DAS Records. *Journal of Geophysical Research: Solid Earth*, 129(1):e2023JB027860, 2024. doi: 10.1029/2023JB027860.

Trabattoni, A., Biagioli, F., Strumia, C., van den Ende, M., Scotto di Uccio, F., Festa, G., Rivet, D., Sladen, A., Ampuero, J. P., Métaxian, J.-P., and Stutzmann, É. From Strain to Displacement: Using Deformation to Enhance Distributed Acoustic Sensing Applications. *Geophysical Journal International*, 235(3): 2372–2384, Dec. 2023. doi: 10.1093/gji/ggad365.

Trabattoni, A., Vernet, C., van den Ende, M., Baillet, M., Potin, B., and Rivet, D. Sediment Corrections for Distributed Acoustic Sensing, Mar. 2024. doi: 10.31223/X5PX19.

Universidad de Chile, t. Red Sismologica Nacional. *International Federation of Digital Seismograph Networks*, 2012. doi: 10.7914/SN/C1.

Viens, L., Bonilla, L. F., Spica, Z. J., Nishida, K., Yamada, T., and Shinohara, M. Nonlinear Earthquake Response of Marine Sediments With Distributed Acoustic Sensing. *Geophysical Research Letters*, 49(21):e2022GL100122, 2022. doi: 10.1029/2022GL100122.

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, I., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., and Contributors, S. . SciPy 1.0–Fundamental Algorithms for Scientific Computing in Python. *arXiv:1907.10121 [physics]*, July 2019.

Yin, J., Soto, M. A., Ramírez, J., Kamalov, V., Zhu, W., Husker, A., and Zhan, Z. Real-Data Testing of Distributed Acoustic Sensing for Offshore Earthquake Early Warning. *The Seismic Record*, 3

(4):269–277, Oct. 2023a. doi: 10.1785/0320230018.

Yin, J., Zhu, W., Li, J., Biondi, E., Miao, Y., Spica, Z. J., Viens, L., Shinohara, M., Ide, S., Mochizuki, K., Husker, A. L., and Zhan, Z. Earthquake Magnitude With DAS: A Transferable Data-Based Scaling Relation. *Geophysical Research Letters*, 50(10): e2023GL103045, 2023b. doi: 10.1029/2023GL103045.

Zhan, Z. Distributed Acoustic Sensing Turns Fiber-Optic Cables into Sensitive Seismic Antennas. *Seismological Research Letters*, 91(1):1–15, Jan. 2020. doi: 10.1785/0220190112.

## Appendix I: gradient-based unwrapping

A somewhat naive approach to phase unwrapping is to detect when the recorded signal exhibits a large discontinuity, which is then interpreted as a phase jump and corrected by adding or subtracting twice the value of the dynamic range. In the context of what was previously discussed at the start of Section 3, this approach amounts to finding the rectangles in Eq. (18). Such an algorithm is currently implemented in the commonly-used NumPy library, and may be a starting point for analysts who wish to restore their saturated DAS recordings. Unfortunately, correctly identifying each phase jump is highly improbable even for modest degrees of saturation, hence this approach is prone to unwrapping errors.

An alternative method circumvents the need for phase jump detection by relying on the continuity of the gradient of the data in the complex plane: a sequence of small increments in the phase angle can be easily tracked, even as it crosses quadrants (e.g. from $+\pi$ to $-\pi$). While the DAS recordings may exhibit discontinuities in the time domain, these are only the result of lifting the phase measurement out of the complex plane (the arg operation in Eq. (4)). Within the complex plane, the signal is continuous, and so it can be reconstructed by estimating the gradient of the signal in the complex plane, converting this into increments of $\hat{\epsilon}$, and integrating (summing) these increments to obtain a signal that is not restricted by the dynamic range.

A first-order algorithm that implements this notion is as follows:

**Algorithm 1 Gradient-based unwrapping**

**Require:** time-series $x$, dynamic range $R$

1: $z \leftarrow \exp(j\pi x R^{-1})$ ▷ Lift $x$ into the complex plane

2: $g \leftarrow \text{diff}(z)$ ▷ Phase difference; see Eq. (2)

3: $\hat{x} \leftarrow \frac{R}{\pi} \int \arg\{g\} dt$ ▷ Reconstruction by integration

4: **return** unwrapped time-series $\hat{x}$

The critical step in this algorithm is the estimation of the gradient of $z$. As in Section 2.1, one can conveniently estimate the gradient in the complex plane through a first-order finite difference operation, but more advanced gradient estimators (such as Kalman filters) can be employed to achieve higher accuracy. Likewise, the integration of $\arg\{g\}$ could take the form of a cumulative sum (if $g$ is expressed as a phase difference) or a higher-order integration scheme (if $g$ is expressed as a phase gradient).
<figure><img src="image_8.png" /><figcaption>**Figure A1** Performance of the gradient-based unwrapping algorithm. The reference signal is scaled by its maximum value, and subsequently multiplied by a constant factor ranging from 1 up to 50. Synthetic recordings are rendered by wrapping the scaled reference signal, which is fed to Algorithm 1. Unwrapping errors start to appear when the maximum amplitude exceeds the dynamic range by a factor 20.</figcaption></figure>

Assuming perfect (infinite-order) differentiation and integration, we can derive a criterion for the onset of unwrapping errors of Algorithm 1. Taking again the example of a monochromatic oscillator, $x(t) = A \cos(2\pi ft)$, and realising that $|\arg \{g\}| < \pi$, we obtain:

$$ |\arg \{g\}| \approx 2\pi^2 AR^{-1} f |\sin(2\pi ft)| \Delta t < \pi \qquad (22) \\ \Leftrightarrow A_{\text{crit}} = \frac{R}{2\pi f \Delta t} $$

The factor $\pi R^{-1}$ is introduced in the first step to acknowledge the conversion from $x$ to $z$. In practice, the differentiation and integration will be imperfect, and so these finite-precision schemes will introduce a proportionality constant, i.e. $A_{\text{crit}} = \beta R (2\pi f \Delta t)^{-1}$, with $0 < \beta < 1$.

As an example, we test Algorithm 1 on a reference signal (recordings of ocean gravity waves) with a characteristic frequency of $f = 0.1$ Hz and a time sampling rate of $(\Delta t)^{-1} = 62.5$ Hz – see Fig. A1. Hence, the theoretical maximum $A_{\text{crit}}R^{-1} \approx 100$. The signal was scaled by its maximum amplitude and multiplied by a factor ranging from 1 to 50 times the dynamic range (i.e., $AR^{-1} \in \{1, 2, 5, 10, 20, 50\}$), followed by wrapping around $R$. After applying Algorithm 1 on the wrapped data, we find

that unwrapping errors start to become prominent after $AR^{-1} > 20$, suggesting $\beta \approx 0.2$. As aforementioned, higher-order differentiation and integration will likely increase $\beta$, possibly approaching $\beta \approx 1$.

The article *An analysis of the dynamic range of Distributed Acoustic Sensing for Earthquake Early Warning* © 2025 by Martijn van den Ende is licensed under **CC BY 4.0**.