# Non-Linear Seismic Velocity Variations Observed During a Seismic Swarm in the Alto Tiberina Low Angle Normal Fault From Ambient Noise Correlation Measurements

Nour Mikhael, Piero Poli, Stéphane Garambois

► To cite this version:

Nour Mikhael, Piero Poli, Stéphane Garambois. Non-Linear Seismic Velocity Variations Observed During a Seismic Swarm in the Alto Tiberina Low Angle Normal Fault From Ambient Noise Correlation Measurements. Journal of Geophysical Research: Solid Earth, 2024, 129 (2), pp.73-100. 10.1029/2023JB028232. hal-04960658

HAL Id: hal-04960658

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

Submitted on 21 Mar 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>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2023JB028232

**Key Points:**

* Velocity variations from seismic coda wave interferometry

* Separation of thermoelastic and hydrological stress induced velocity variation from tectonic changes

* Seismic swarm triggered by aseismic slip favored by the presence of fluids

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

P. Poli and N. Mikhael,
piero.poli@unipd.it;
nour.mikhael@univ-grenoble-alpes.fr

**Citation:**

Mikhael, N., Poli, P., & Garambois, S. (2024). Non-linear seismic velocity variations observed during a seismic swarm in the Alto Tiberina low angle normal fault from ambient noise correlation measurements. *Journal of Geophysical Research: Solid Earth*, 129, e2023JB028232. https://doi.org/10.1029/2023JB028232

Received 9 NOV 2023
Accepted 6 FEB 2024

**Author Contributions:**

**Conceptualization:** Piero Poli,
Stéphane Garambois

**Formal analysis:** Nour Mikhael,
Piero Poli, Stéphane Garambois

**Funding acquisition:** Piero Poli

**Investigation:** Nour Mikhael

**Methodology:** Nour Mikhael, Piero Poli,
Stéphane Garambois

**Project administration:** Piero Poli

**Resources:** Piero Poli

**Software:** Nour Mikhael

**Supervision:** Piero Poli,
Stéphane Garambois

**Validation:** Piero Poli,
Stéphane Garambois

**Visualization:** Nour Mikhael

**Writing – original draft:** Nour Mikhael,
Piero Poli, Stéphane Garambois

© 2024. The Authors.

This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

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

## Non-Linear Seismic Velocity Variations Observed During a Seismic Swarm in the Alto Tiberina Low Angle Normal Fault From Ambient Noise Correlation Measurements

Nour Mikhael<fnref n="1" />,*, Piero Poli<fnref n="1" />,<fnref n="2" />, and Stéphane Garambois<fnref n="1" />

**Abstract** From seismic interferometry, we investigate the strain sensitivity to seismic velocity variations related to a seismic swarm activity that occurred in 2013 along the Alto Tiberina low angle normal fault. We compute daily auto-correlation functions of ambient noise recorded at seismic stations located in the vicinity of the fault over the course of 10 years. Using the stretching technique, we compute daily velocity variations smoothed over a period of 100 days with a time lapse approach. Through the application of an optimization procedure based on synthetic modeling, we separate the non-tectonic, thermoelastic and rain induced velocity variations, from the tectonic components. Consequently, we unravel a significant velocity drop of 0.033% coinciding with the swarm occurring at seismogenic depth (3–5 km). Additionally, the time evolution of the velocity changes shows a direct relationship with the strain rate rather than the strain indicating a non-linear behavior of the crust induced by aseismic slip. The deduced strain sensitivity, exhibiting an order of magnitude comparable to that observed within volcanic settings, confirms this non-linear behavior and suggests the presence of pressurized fluids at depth.

**Plain Language Summary** The driving mechanism of seismic swarms is still unclear. Seismic swarms are bursts of small earthquakes without a clear triggering mainshock. Understanding the physical processes driving these sequences is crucial for earthquake hazard assessment. To that scope, we focus in this study on a small seismic swarm that occurred in 2013 in northern Apennines, Italy. We monitor the properties of crustal rocks using continuous recordings of ambient seismic noise recorded at seismic stations located in the vicinity of the Alto Tiberina low angle normal fault over the course of 10 years. With this approach, we demonstrate the ability to retrieve quantitative information about the fault's state and to indirectly infer that the triggering of the swarm is mainly related to the presence of fluids in the crust.

### 1. Introduction

Earthquake swarms are sequences of relatively small to moderate earthquakes lasting from hours to months without a clear mainshock driving them, conversely to classic mainshock-aftershock sequences (Hill, 1977; Lohman & McGuire, 2007; Passarelli et al., 2018). Swarms have been observed in different tectonic settings and generally result from some external forcing which drives the occurrence of seismicity (De Barros et al., 2020; Lohman & McGuire, 2007; Passarelli et al., 2018). In volcanic and geothermal regions, the triggering of the seismic swarm is often associated with magma intrusion and fluid migration (Bourouis & Bernard, 2007; Poli et al., 2022). These dynamic and transient processes increase pore pressure, which reduces the strength of the fault and favors the occurrence of earthquakes. Along transform faults boundaries, it is suggested that the driving mechanism that triggers seismic swarms is aseismic creep (Lohman & McGuire, 2007; Roland & McGuire, 2009). More recently, a combination of fluids and aseismic slip has also been proposed from observations of induced earthquake activity (Guglielmi et al., 2015) and of natural swarms (De Barros et al., 2020).

The dynamics of microseismicity (Bourouis & Bernard, 2007; De Barros et al., 2020) and geodetic observations (Cervelli et al., 2001; Gualandi et al., 2017) are the main source of information about the physical processes involved during swarms. The spatiotemporal evolution of seismicity is often used to distinguish in between linear and diffusive evolution of seismicity, behaviors associated with aseismic deformation and fluid migration respectively (Fischer & Hainzl, 2021; Kato et al., 2012; S. A. Shapiro et al., 1997). However complex seismic behaviors, documented in both in-situ controlled experiments (De Barros et al., 2018; Guglielmi et al., 2015) and from natural swarms analysis (De Barros et al., 2020), suggest that there is an interplay of multiple driving

<footnote n="1">¹University Grenoble Alpes, University Savoie Mont Blanc, CNRS, IRD, University Gustave Eiffel, Grenoble, France,</footnote>

<footnote n="2">²Dipartimento di Geoscienze, Università degli Studi di Padova, Padova, Italy</footnote>
Writing – review & editing:

Nour Mikhael, Piero Poli,
Stéphane Garambois

processes at work indicating that aseismic slip is potentially aided by the migration of deep pressurized fluids. The enhancement of aseismic slip through fluid migration can occur due to a fault weakening attributed to an increase in the pore pressure and to the creation of new fractures. This highlights the complexity of adopting the spatiotemporal evolution of seismicity as proxy for determining the triggering mechanisms of swarms. Resolving the underlying causes behind earthquake triggering often requires the integration of multiple geophysical measurements. Thus, emphasizing the importance of conducting complementary analyses that, when integrated with seismicity migration studies or others, can provide a more comprehensive insight into the driving mechanisms behind seismic activity. More difficult is the direct resolution of aseismic slip with geodetic data only (Gualandi et al., 2017; Lohman & McGuire, 2007; Roland & McGuire, 2009), due to the lack of resolution of geodesy to monitor earthquake sequences with cumulative magnitude lower than ~5.

Little is known about quantitative physical conditions of rocks favoring the occurrence of swarms rather than more classic mainshock-aftershocks sequences. General models describing the fault systems promoting the occurrence of swarms exist and suggest that both fluids and complex crack distribution play a role in the generation of these peculiar seismic sequences (Bachura et al., 2021; Hill, 1977; Passarelli et al., 2018; Vidale & Shearer, 2006).

A whole body of experimental setups are nowadays available to assess rocks characteristics including the dynamic moduli such as shear, bulk, and Young's modulus, the mechanical hysteresis, the resonance and so forth (Ostrovsky & Johnson, 2001; Riviere et al., 2014). By assessing the changes in the elastic modulus (e.g., variation of seismic velocity) of rocks that are subjected to stress perturbations, the non-linear stress sensitivity of rocks can be highlighted and permits to quantify the inherent physical conditions (elastic properties, damage, confining pressure) of rock samples (Johnson & Jia, 2005; Ostrovsky & Johnson, 2001; Riviere et al., 2014). It has been shown, for example, that complex material shows strong non-linearity and large modulus variations under small stress perturbation (Johnson & Jia, 2005; Riviere et al., 2014). This sensitivity to stress perturbation is further enhanced by the reduction of confining pressure, or by the increment of material damage during fatigue experiments (Riviere et al., 2014). Changes in the pore pressure within the crust significantly alter the sensitivity to stress perturbation: an increase in the pore pressure tends to reduce the effective confining pressure rendering rocks more susceptible to failure and fracturing. Such reduction also facilitates fluid migration through the enhanced fracture permeability thus intensifying the crack density within the crust and subsequently increasing the sensitivity to stress perturbation (Taira et al., 2018).

Earth rocks are known to respond non-linearly to strain perturbations (Agnew, 1981; Yamamura et al., 2003). Changes in their physical properties can be assessed through the monitoring of seismic velocity variation ($dv/v$), however the continuous monitoring of seismic $dv/v$ mostly relies on active source experiments (Yamamura et al., 2003), which makes the continuous monitoring of faults systems prohibitive or limited to short periods. Ambient seismic noise correlation is nowadays routinely used to monitor seismic velocity changes in various environments and at various scales. However, the interpretation of velocity changes measured with ambient seismic noise correlation in real fault systems remains challenging, mainly because tectonic induced changes are small, unless after the occurrence of large magnitude earthquakes (e.g., $M_w$ 6+) (Brenguier, Campillo, et al., 2008; Brenguier et al., 2014; Poli et al., 2020). Moreover, these minor changes are often hidden by velocity variations induced by meteorological processes such as rain or temperature changes (Barajas et al., 2021; Poli et al., 2020; Rivet et al., 2015; Wang et al., 2017).

A recent study conducted by Sheng et al. (2022) employed velocity variation measurements to discern velocity changes associated with a previously undetected slow slip event at depth. This particular slow slip event falls below the detection threshold attainable through geodetic means and yet was inferred using velocity measurements. Furthermore, the work of Rivet et al. (2011) demonstrated the feasibility of detecting slow slip events through $dv/v$ monitoring establishing a precedent for such $dv/v$ applications. Therefore, in addition to geodetic measurements and seismicity evolution analysis, velocity measurements serve as a complementary source of information shedding light onto the underlying physical processes governing seismic swarms and their potential distinction.

Here, we investigate seismic velocity variations associated with a seismic swarm that took place in the northern Apennines (Figure 1a) (Essing & Poli, 2022; Gualandi et al., 2017). Historically, this tectonically active region witnessed the occurrence of large earthquakes with magnitudes higher than $M_w$ 5 (De Gori et al., 2015; Gualandi et al., 2017; Vadacca et al., 2016). Continuous microseismicity ($M_l$ < 2.3) has also been reported on the Alto
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Elevation map of the northern Apennines (see inset) showing the spatial distribution of the seismic swarm of 2013–2014 (colored dots). Black triangles refer to the seismic stations used in this study. Blue inverted triangles refer to GPS stations. The main faults, the ATF (east) and the GF (central), are mapped in red. The star indicates the location of largest event of the swarm along with a beachball to represent its focal mechanism. (b) Relative distance evolution between the two GPS stations MVAL and ATBU localized on opposite compartments of the GF fault, in gray, and daily rate of seismic events in red.</figcaption></figure>

Tiberina fault (ATF) along with swarm like activity along the Gubbio fault (GF) (De Gori et al., 2015; Essing & Poli, 2022; Gualandi et al., 2017; Vadacca et al., 2016). Our focus is the analysis of the seismic swarm activity of 2013 surrounding the GF, which was accompanied by a small surface deformation expressed by a fault slip (Figure 1b). The GPS data used to estimate this surface deformation can be found at Blewitt et al. (2018). The activity started in 2013 and lasted until the end of 2015. The largest event in this sequence is of a magnitude $M_w$ 3.9 (Figure 1a) (Essing & Poli, 2022), while the total moment released by the swarm is equivalent to a magnitude $M_w$ 4.3 earthquake (Gualandi et al., 2017). Through the exploit of temporal $dv/v$ we resolve a 0.033% velocity drop associated with the seismic swarm and its related geodetic deformation (Figure 1b). This is achieved by deconvolving environmental processes that are simultaneously contributing to the $dv/v$ (Figure 2), namely thermoelastic stress and loading related to climatic forcing. Our observation suggests that the presence of pressurized fluids at depth induced aseismic slip which in return triggered the seismic activity in the ATF system.

## 2. Data and Methods

Temporal seismic velocity variations that occur in the crust can be detected through coda wave interferometry (Sens-Schönfelder & Wegler, 2006; Wegler & Sens-Schönfelder, 2007). This method is based on the reconstruction of the empirical Green's function (EGF) between seismic stations through ambient noise correlation (Campillo & Paul, 2003; N. M. Shapiro & Campillo, 2004; N. M. Shapiro et al., 2005). The EGF describes the Earth's impulse response between a pair of seismic stations. Its reconstruction is contingent on the convergence of the correlations as well as the temporal stability and the homogeneous spatial distribution of the ambient noise sources. Though, it has been proven, that for monitoring purposes, only the condition on the temporal stability of the noise sources is required (Hadziioannou et al., 2009).

Three distinct correlation techniques, namely the crosscorrelation (CC), the autocorrelation (AC) and the single station crosscorrelation (SC), yield an estimation of the EGF (Hobiger et al., 2014, 2016). Within the literature, it has been documented that the AC method has been effectively used for retrieving $dv/v$ associated with large earthquakes (Brenguier, Shapiro, et al., 2008), for monitoring annual and multi-annual seasonal variations (Li & Ben-Zion, 2023; Poli et al., 2020) leading up to detecting small $dv/v$ related to seismic swarms (Maeda et al., 2010). Furthermore, the AC method facilitates the precise assignment of $dv/v$ to individual seismic stations avoiding the complexities associated with the spatial allocation of $dv/v$ in the context of the CC method, making it the technique of choice for this study. A comparative analysis of these three methods conducted by Hobiger et al. (2012) demonstrates that the results are frequency dependent and show that better outcomes are observed from the AC approach for frequencies exceeding 0.5 Hz. This result leads to focus our computations on a
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Stack of relative velocity variations measured from auto-correlation of seismic noise recorded at several stations in the frequency band [0.5–1] Hz for several coda lapse times along with their associated standard deviation in light orange. The red shaded rectangle highlights the time period of the seismic swarm.</figcaption></figure>

frequency range spanning from 0.5 to 1 Hz for the *dv/v* investigation. This is depicted in Figure S1 in Supporting Information S1 where we compare *dv/v* computed within two frequency bands ([0.1–0.5] Hz and [0.5–1] Hz) and observe a clear enhancement in both the signal-to-noise ratio and the stability of the *dv/v* within the highest band. It's worth noting that restricting the analysis to a single higher frequency band does not necessarily limit the depth extension of the *dv/v* analysis when employing a time-lapse approach. Further discussion on the matter is subsequently presented.

## 2.1. Auto-Correlation Functions

We use continuous seismic data recorded at five stations located in the vicinity of the Gubbio fault (Figure 1a) belonging to the INGV (Istituto Nazionale di Geofisica e Vulcanologia) network (INGV, 2005). We select 10 years of data, from the 1st of January 2011 to the 31st of December 2020, comprising the occurrence period of the seismic swarm (Gualandi et al., 2017). The processing of the records and the retrieval of the EGFs are similar to the processing described in Poli et al. (2020). We compute the three components (EE, NN, ZZ) auto-correlation functions (ACFs) on a daily basis for all stations. For that, we divide one-day traces into ten-min segments with a five-min overlap. After removing the instrumental response, we apply a bandpass filter between 0.5 and 1 Hz. To remove spurious spikes and contributions from large earthquakes, we remove samples having an energy exceeding the daily mean of the energy plus three standard deviations. Low magnitude earthquakes at such frequencies display a low signal to noise ratio, and thus do not contaminate the filtered signal. At this point, if our daily trace does not contain more than 20 hr of data available, the trace is rejected from the data set. To further reduce transient signals, we carried out a one-bit normalization (Poli et al., 2020). Then, we compute ACFs for all ten-min segments and average them to obtain a single station daily ACF (Figure S2 in Supporting Information S1). A reference ACF representative of the

background crustal state is computed for each station by averaging all daily ACFs of those stations (Figure S2 in Supporting Information S1). To improve the temporal stability of the ambient noise sources, and consequently the *dv/v* measurements, we stack the daily ACFs over 100 days with a 99-day overlap.

## 2.2. Stretching

The velocity changes are then measured from the stacked ACFs with the stretching technique (Poupinet et al., 1984; Snieder et al., 2002). This method treats the daily ACFs as a linearly stretched or a compressed version of a reference ACF. In practice, the stretching is done by resampling the ACF with a sampling rate increased by a factor one minus the stretching coefficient. Assuming a homogeneous velocity perturbation, this stretching coefficient equals a velocity change (Lobkis & Weaver, 2003; Zhao et al., 2010). The similarity between the reference and current ACF is evaluated through their correlation coefficient for a given stretching coefficient. The stretching coefficient that maximizes the correlation coefficient is the one that best fits the data and represents the seismic velocity change (Sens-Schönfelder & Wegler, 2006, 2011).

For all of the stations, each daily ACF was stretched and compressed using a regular grid search with 2000 trials between -0.8% and 0.8%. If the *dv/v* value does not have a correlation coefficient higher than 0.75, the value is excluded from the analysis. We then average the *dv/v* over the three components using a weighted approach based on their squared correlation coefficients obtained from stretching.

The *dv/v* are measured in a series of time windows across the coda waves (Figure 2), in order to scan different depths of the fault system (Poli et al., 2020). Numerical simulations (Obermann et al., 2013) and real data investigations (Obermann et al., 2014) illustrate that deeper portions of the crust are sampled by later coda time windows. We investigate five coda windows of 20s length to avoid cycle skipping considering a frequency band
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Daily maximum temperature. (b) Daily rainfall. Velocity variations in black for coda lapse time 10–30s along with the synthetic velocity variations according to the temperature model in orange (c) and the rain model in blue (e). Velocity variations in black for coda lapse time 34–54s along with the synthetic velocity variations according to the temperature model in orange (d) and the rain model in blue (f). (g) Correlation coefficient between synthetic velocity variations and measured velocity variations as a function of coda lapse time window. Blue and orange dots refer to the synthetic cases according to the consideration of rain and temperature models respectively, in an independent approach.</figcaption></figure>

of [0.5–1] Hz within lag times ranging from 10 to 54s. The choice of the coda time windows is guided by previous studies (Obermann et al., 2013; Poli et al., 2020) and signal-to-noise ratio considerations.

## 3. Velocity Variations

The relative velocity changes obtained from the analysis of the ambient noise auto-correlation are spatially averaged over the stations (Figure 2). Several temporal patterns can be distinguished. We observe dominant yearly seasonal variations for the first coda lapse time 10–30s, while this pattern vanishes in the later coda windows. A similar behavior has been previously observed in the Apennines (Poli et al., 2020) and suggests a dominant effect of thermoelastic stress in modulating the *dv/v* for shallow crustal layers (Richter et al., 2014). On the other hand, the later coda, which is sensitive to deeper layers, show smaller velocity variations, with longer periods, likely related to long term hydrological effects (Poli et al., 2020). No direct evidence of seismic velocity changes is observed during the swarm activity (Figures 1b and 2).

### 3.1. Temperature-Induced *dv/v* Model

To improve our understanding of the effect of the thermoelastic stress in controlling the observed *dv/v*, we generate temperature-induced synthetic *dv/v* time series. To achieve this, we procure daily temperature data sets from the NASA—LaRC—POWER project (https://power.larc.nasa.gov/) data catalog found at Stackhouse (2021). The data set has a spatial resolution of 55 km × 77 km, encompassing the entire spatial distribution of the seismic stations. In this context, the use of such a broad resolution is justified considering that the retrieval of synthetic data is intended for station-averaged *dv/v*. Then, a moving average of 100 days is applied to the recorded daily temperature (Figure 3a) in accordance to the temporal stacking of the ACFs. Subsequently, we implement the transfer Equation 1 to our station-averaged *dv/v* measurements leading to the derivation of the temperature-induced synthetic *dv/v* time series, as documented in previous studies (Barajas et al., 2021; Rivet et al., 2015). In Equation 1, < > denotes the average over time, *T* the 100-day smoothed temperature variation, *cov* the covariance function and *var* the variance function.
$$ \frac{dv}{v_{\text{syn}}}(t) = \left\langle \frac{dv}{v}(t) \right\rangle + \left( \frac{\operatorname{cov}(\frac{dv}{v}(t), T(t))}{\operatorname{var}(T(t))} \right) \cdot T(t) \quad (1) $$

The results show that for early coda, velocity changes are controlled by temperature (Figure 3c). For late coda, temperature-based models cannot explain the observed velocity changes (Figure 3d). This confirms that temperature plays a primary role in controlling the velocity changes at shallow depths (Poli et al., 2020; Richter et al., 2014).

## 3.2. Rain-Induced dv/v Models

To explore if the tectonic stress induced variations are hidden by non-tectonic contributions, and particularly by hydrological effects, we use the synthetic dv/v model proposed by Barajas et al. (2021). For that, we computed rain-induced dv/v from estimates of the accumulation of water level in the crust with respect to time, derived from rainfall data. Again, since this approach is based on a station-averaged analysis we produce daily rainfall data sets from the NASA—LaRC—POWER project as we did for the temperature.

This synthetic dv/v retrieval approach assumes a recharge-discharge ongoing process, where the daily rainfall infiltrates the crust and contributes to the groundwater recharge. As for the discharge process, we assume that drainage occurs in terms of a Torricelli efflux velocity (Fiorillo, 2011). It should be noted that such model does not consider external factors that may vary the discharge or/and the recharge of the reservoir such as evapotranspiration. Thus, the discharge rate can be expressed in terms of the drainage velocity as stated by Torricelli's law:

$$ \frac{dQ}{dt} = A_s \sqrt{2gh} + R, \quad (2) $$

where Q is the volume of water, g is the gravity, h is the water level in the crust, R is the external source supplying the reservoir with water and A_s is the area through which the water is discharged. By expressing Q as the product of the area of the bottom A_B with the height of the water column h we can re-write Equation 2 as:

$$ \frac{dh(t)}{dt} = k \sqrt{h(t)} + r, \quad (3) $$

with k = A_s/A_B√2g and r = R/A_B. If we discretize Equation 3 with a first order approximation we end up with the following expression:

$$ h_{i+1} = h_i + (t_{i+1} - t_i) \cdot k \sqrt{h_i} + (t_{i+1} - t_i) \cdot r_i. \quad (4) $$

Since our time resolution is of one day, we can reduce Equation 4 to:

$$ h_{i+1} = h_i + k \sqrt{h_i} + r_i. \quad (5) $$

The synthetic model of velocity variations can be estimated as done for the temperature using Equation 1 by substituting the temperature parameter T with the water level estimates h. If the observed dv/v are driven by rainfall, there exists a value of the parameter k that optimizes the misfit between observed and synthetic dv/v. In order to find this optimal k value, we compute synthetic dv/v models for various k values ranging from 0 to 0.5 with a sampling step of 2 · 10⁻⁴. We later assess the model fit to the measured data with the misfit σ(k)² (6). The best fit corresponds to the k value that gives the minimal misfit.

$$ \sigma(k)^2 = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{dv}{v}(i) - \frac{dv}{v_{\text{syn}}}(i, k) \right)^2. \quad (6) $$
<table><caption>**Table 1**
σ(k)² Values Between Synthetic and Observed dv/v for All Coda Lapse Time Windows and Each of the Three Models Used</caption><thead><tr><th>Coda windows</th><th>10–30s</th><th>16–36s</th><th>22–42s</th><th>28–48s</th><th>34–54s</th></tr></thead><tbody><tr><td>Model 1</td><td>5.11 · 10<sup>-4</sup></td><td>2.27 · 10<sup>-4</sup></td><td>2.48 · 10<sup>-4</sup></td><td>2.26 · 10<sup>-4</sup></td><td>2.75 · 10<sup>-4</sup></td></tr><tr><td>Model 2</td><td>12.52 · 10<sup>-4</sup></td><td>4.38 · 10<sup>-4</sup></td><td>2.93 · 10<sup>-4</sup></td><td>2.22 · 10<sup>-4</sup></td><td>2.32 · 10<sup>-4</sup></td></tr><tr><td>Model 3</td><td>4.82 · 10<sup>-4</sup></td><td>1.83 · 10<sup>-4</sup></td><td>1.88 · 10<sup>-4</sup></td><td>1.34 · 10<sup>-4</sup></td><td>1.86 · 10<sup>-4</sup></td></tr><tr><td>Model 4</td><td>4.56 · 10<sup>-4</sup></td><td>1.92 · 10<sup>-4</sup></td><td>1.43 · 10<sup>-4</sup></td><td>1.38 · 10<sup>-4</sup></td><td>1.70 · 10<sup>-4</sup></td></tr></tbody></table>

Note. Models 1, 2, 3, and 4 refer to temperature based, rain based, combined model on averaged dv/v and combined models on single station dv/v respectively.

The optimized models show that for early coda, rain-based models cannot account for the large annual fluctuations visible in the *dv/v* (Figure 3e) while for late coda, an enhanced performance is observed (Figure 3f). This further confirms that shallow layers are affected by thermal stress induced variations and deep layers are more likely to be sensitive to rain induced stress changes (Poli et al., 2020).

This trend is quantified using the correlation coefficient (cc) evolution computed between the synthetic and the measured *dv/v* variations with respect to the coda lapse time window (Figure 3g). We observe a high cc value of 0.78 between the temperature model and the measured *dv/v* for the first coda window that drastically drops afterward to reach a final cc value of 0.31. In the meantime, cc values between the rain models and the measured *dv/v*, despite the low cc values exhibited, a significant increase from 0.18 for the first coda lapse time to a cc value of 0.46 for the last coda lapse times is noticeable (Figure 3g). The inadequacy of both modeling approaches in fully explaining the observed *dv/v* for the late coda implies the existence of either an unaccounted for third component driving the *dv/v* or a complex combination of multiple processes governing the velocity changes.

### 3.3. Tectonic Velocity Variation

To take into account the simultaneous coda dependent contributions of both temperature and rainfall, we consider a generalized model depicting a linear relationship between the measured *dv/v* and the combined effects, as in Wang et al. (2017):

$$ \frac{dv}{v_{\text{syn}}} (t) = a \cdot T(t) + b \cdot h(t) + c \cdot t, \qquad (7) $$

where *a* and *b* are the controlling coefficients of the temperature and water level height respectively, and *c* is the coefficient for the linear trend if any exist.

In light of the well-documented variation in precipitation patterns within mountainous regions due to orographic influences (Daly et al., 1994), we implement this generalized model in two distinct manners. We designate them as “Model 3” and “Model 4” respectively, while bearing in mind that “Model 1” and “Model 2” correspond to the temperature-induced synthetic *dv/v* and the rain-induced synthetic *dv/v*, respectively presented earlier. Model 3 considers regional climatological factors while Model 4 addresses the potential climatological variations that may manifest at the level of individual stations. To that end, we use a higher resolution environmental data set, for both the daily temperature and the daily rainfall, from the Joint Research Centre Data Catalog—EMO (Thiemig et al., 2022). This data set is characterized by an enhanced spatial resolution of 1.5 km × 1.5 km.

To provide further clarification, in Model 3, we apply Equation 7 to the regional environmental data set and the station-averaged *dv/v*, resulting in a single generalized synthetic *dv/v*. Conversely, in Model 4, we independently apply Equation 7 to *dv/v* obtained per station. This approach mitigates the influence of local climatological variabilities. Subsequently, we obtain several generalized synthetic *dv/v*, one for each station, which are then stacked spatially.

To assess the precision of our models, we conduct a comparative analysis by measuring the misfit (Equation 6) between the observed and synthetic *dv/v* (Table 1). Our findings reveal that a generalized model, considering the simultaneous impact of both precipitation and temperature (Model 3 and Model 4), outperforms individual models (Model 1 and Model 2). Furthermore, when contrasting the station-averaged approach (Model 3) with the
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Relative velocity variations measurements (black) and optimized synthetic velocity variations models (green) in the frequency range [0.5–1] Hz for several coda lapse times.</figcaption></figure>

station-wise approach (Model 4), we observe a high similarity within the misfit values. Consequently, we opt to adopt Model 4, given its suitability for station-wise analyses when required (Figure 4). The *k* parameter that yielded said optimized misfits is shown in Figure S3 in Supporting Information S1. Taking into account the negligible influence of precipitation-induced perturbations on *dv/v* within early time lags, as previously indicated, it is advisable to exclusively interpret *k* within late time lags. We observe a consistent alignment of the *k* parameter among seismic stations during late coda windows, with no large spatial deviations.

When concurrently accounting for the influence of temperature and rainfall, it becomes evident that the synthetic *dv/v* closely align with the observed *dv/v* across all coda windows (Figure 4). Interestingly, for the late coda window (Figure 4e), we observe notable discrepancies between the synthetic *dv/v* and the observed *dv/v*. Particularly, a velocity drop is observed around the year 2014 coinciding with the swarm activity (De Gori et al., 2015; Essing & Poli, 2022; Gualandi et al., 2017), wemains unexplained by the synthetic *dv/v*. To better focus on this velocity change, we removed the non-tectonic component, estimated from our generalized synthetic approach, in order to focus on the residual tectonic sensitivity of *dv/v* measurements (Figure 5). This process was again performed on a single station basis, where the residuals from each station were first retrieved and then averaged providing a general tectonic residual.

To evaluate the stability of our residual *dv/v* (Figure 5a) and estimate a confidence interval, we use the jackknife resampling method (Figure S4 in Supporting Information S1) and assess the number of stations contributing to the averaging process. The jackknife resampling involves the systematic re-evaluation of a selected statistic across a collection of observations while, at each iteration, randomly excluding specific observations from that ensemble. In our case, the statistic of interest is the mean, the ensemble of observations is the set of residual *dv/v* per stations, and the total number of iterations is set to 1000.

We notice that during most iterations, we consistently obtain a stable average residual *dv/v*, exhibiting a high correlation coefficient (>0.7) with the target residual *dv/v* (Figure S4 in Supporting Information S1). Nevertheless, during two distinct time instances, namely early 2018 and late 2019, we observe a considerable increase in
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Residual velocity variations (black) along with their standard deviation (light orange) in the frequency range [0.5–1] Hz for the coda window 34–54s. The shaded rectangles correspond to time periods exhibiting large velocity drops. (b) Number of missing stations from the *dv/v* averaging process in gray dots.</figcaption></figure>

the standard deviation compared to other time periods (Figures S4 and S5a in Supporting Information S1) which implies that the interpretation of *dv/v* within those particular time periods should be done with caution due to high uncertainty.

When examining deviations within the synthetic *dv/v* with respect to the observed *dv/v* (Figures 4e and 5a), we identify three main time periods characterized by a sharp velocity drop spanning the duration of few weeks in around 2015, 2018, and 2019 (indicated in gray). We infer that those velocity drops observed within the gray-highlighted periods exhibiting high uncertainties arise from a decrease in the count of contributing stations to the averaging process (Figure 5b). Additionally, a long-term gradual decrease of velocity is observed from 2013 to 2014. This velocity drop observed during the red-highlighted period coinciding with the occurrence of the seismic swarm, exhibits comparatively low uncertainties and does not correlate with station coverage gaps (Figure 5).

Consequently, this residual *dv/v* in the late coda window (Figure 5a) is mainly characterized by a decrease in *dv/v* that reaches its maximum end of 2013 beginning of 2014 concurring with the time period of the seismic swarm. This yields a total velocity drop of 0.033% ± 0.009%, which is then followed by almost a 1-year recovery to the original state.

## 4. Discussion

Our detailed analysis of seismic velocity variations measured nearby the ATF fault system, revealed a dominant contamination of *dv/v* time series by environmental processes, such as rain and temperature variations (Figure 3). The modeling of *dv/v* induced by these external processes highlights how early coda waves are mainly affected by shallow temperature variations, inducing thermal stress evolution (Poli et al., 2020; Richter et al., 2014). The *dv/v* measured from later coda waves appear more sensitive to long periodic variations of rainfall, in agreement with previous results (Poli et al., 2020).

The derived synthetic *dv/v* models (Figure 4) permitted to obtain residual time series, which show an excess of negative velocity variations for the late coda windows, not explained by environmental factors (Figure 5), but rather correlated with aseismic slip and increase of microseismicity observed in the period 2013–2014. The largest tectonic induced velocity variations are only observed in the later coda (Figures 4e and 5).

Accurate depth characterization of the observed velocity variation requires coda wave sensitivity kernels obtained through an inversion process (Barajas et al., 2022). However, a preliminary estimation of the depth can be assessed by examining the distinct behaviors within each coda window. Early coda lapse times primarily reflect the influence of thermal stress changes (Figure 3). This thermal effect is known to rapidly decay as function of depth and is confined to the shallow layers of the crust (Richter et al., 2014). Furthermore, the *dv/v* residuals within the early coda window exhibit no discernible tectonic signature after removal of the environmental effects. Conversely, a tectonic signature is observed exclusively in the late coda window indicating that the changes were generated within a restricted volume surrounding seismic depth (3–5 km) (Gualandi et al., 2017; Obermann et al., 2013; Poli et al., 2020). The dichotomy of behaviors observed in the coda windows therefore allows for a preliminary qualitative depth characterization. Thus, our measurements of *dv/v* over different coda lapse times, together with detailed modeling of external environmental processes, permitted to resolve velocity variations related to a small (equivalent magnitude of *M*<sub>w</sub> 4.3) seismic swarm.

To incorporate the tectonic signature observed in the residual *dv/v* time series (Figure 5) into our synthetic modeling, we propose a modification to Equation 7. This adjustment involves the introduction of a tectonic term associated with the 2013–2014 seismic swarm based on the deformation rate. The observed velocity reduction shows a strong correlation with the deformation rate rendering it a favorable proxy for generating synthetic *dv/v* time series that are driven by seismic processes. This fifth and last model “Model 5” is as follows:
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Relative velocity variations measurements (black), optimized synthetic velocity variations derived from model 4 (green) and corrected synthetic velocity variations (maroon) in the frequency range [0.5–1] Hz for the coda window 34–54s. The dashed curve represents the distance rate. (b) Residual velocity variations corresponding to each of the synthetic velocity variations for the coda lapse time 34–54s. (c) Relative distance evolution between the two GPS stations MVAL and ATBU localized at opposite sides from the GF in gray and daily rate of seismic events in red. The dashed curve represents the ramp fit of the distance. The gray shaded rectangles in (a, b) delimit a previously accounted for velocity drop.</figcaption></figure>

$$ \frac{dv}{v_{\text{syn}}} (t) = a \cdot T(t) + b \cdot h(t) + c \cdot t + d \cdot D_r(t), \qquad (8) $$

with *d* being the amplitude of the velocity drop associated with the seismic swarm and *D*<sub>r</sub>(*t*) the deformation rate.

The estimation of the deformation rate is accomplished through a threefold procedure. Initially, we approximate the time fluctuation of the relative distance between the two GPS stations using a ramp-shaped function (Figure 6c). Then, we compute the derivative of this fitted ramp function, resulting in a boxcar function. To enhance the cohesion among all constituent elements of the *dv/v* (temperature, rain and seismic deformation), we apply a smoothing to the edges of this boxcar function by convolving it with a Hanning window with a window length set to 100 days. The choice of a 100-day window length is made to account for the original stacking process applied to the ACFs. Consequently, the resulting output corresponds to the target deformation rate (Figure 6a).

The impact of the adjustment made within the modeling process is depicted in Figure 6a for the coda window 34–54s and in Figure S5 in Supporting Information S1 for all coda windows. We conduct a comparative analysis between the synthetic *dv/v* time series generated by Model 4 and Model 5 and the observed *dv/v* time series. We notice that the previously unaccounted for excess velocity decrease around the 2013–2014 period, which was not
captured by Model 4, is now accurately represented by Model 5. Furthermore, we note that the residual *dv/v* obtained from Model 5 no longer exhibit a velocity decrease during the occurrence of the seismic swarm, instead they fluctuate around zero (Figure 6b). This transition from Model 4 to Model 5 also results in a reduction of the misfit from 1.70 · 10⁻⁴ to 1.38 · 10⁻⁴. This highlights the critical importance of incorporating an additional seismic term to enhance the accuracy and reliability of synthetic *dv/v* modeling that faithfully mirrors the true behavior of the measured *dv/v*.

This approach also enables the determination of station-specific velocity drop magnitudes. Our analysis reveals that the velocity drop magnitudes for stations IV.ATFO, IV.MURB, IV.ATVO, IV.ATPC, and IV.ATMI are 0.027%, 0.020%, 0.011%, 0.008%, and 0.035%, respectively with a total average of 0.020%. It is anticipated and reasonable to observe a decrease in the velocity drop magnitude when stations are located farther away from the fault. However, the furthest station, station IV.ATMI, displays an unconventional pattern, having the highest drop magnitude. This anomaly may be attributed to site effects resulting from the location of the IV.ATMI station within a valley (Boschelli et al., 2021). This site effect involves the entrapment of seismic waves in the valley which results in a higher level ground motion consequently yielding a larger velocity reduction.

The time evolution of velocity variation reveals a close relationship with the fault opening observed at surface by geodetic data, with peak *dv/v* occurring at peak deformation rate (Figure 6). This observation is similar to the velocity variation observed during slow deformation experiments conducted in laboratory (Scuderi et al., 2016) and slow slips in subduction zones (Rivet et al., 2011). As in previous studies (Rivet et al., 2011), we notice that this behavior indicates non-linear (or non-elastic) response of crustal material under low strain rates induced by aseismic slip.

The time evolution of *dv/v* also correlates with the occurrence of substantial number of earthquakes from the beginning of 2014 (Figure 6). While velocity changes are primary induced by dynamic strain generated by large earthquakes (Brenguier, Shapiro, et al., 2008), it is unlikely here that the low shaking induced by the low magnitude events composing the swarm activity (Essing & Poli, 2022) is responsible for the observed velocity change (Figure 6). We thus propose that the *dv/v* is induced by aseismic slip.

To assess how the aseismic slip (Gualandi et al., 2017) affects the velocity variation, we model the volumetric strain ε. For that, we used a fault slip model for the 2013–2014 deformation episode based on the work of Gualandi et al. (2017) and the Coulomb software (Toda et al., 2011) to estimate the spatial distribution of the volumetric strain assuming an elastic half space with uniform isotropic elastic properties (following Okada (1992)). We compute the volumetric strain at the surface (Figure S6a in Supporting Information S1), along different depth levels and along a cross section perpendicular to the fault (Figure S6b in Supporting Information S1). This strain model is a basic representation founded on a linear elastic approach. To achieve a precise 3D estimation of the strain, the incorporation of kernel sensitivity models is required which is outside of the scope of this study. While the linear elastic strain model has proven effective for large earthquakes, note that it is applied here to a seismic swarm.

All stations show a drop in *dv/v*, whether they are located in the dilatation or compression areas predicted by the linear elastic deformation model (Figure S7 in Supporting Information S1). This lack of correlation suggests the presence of damage due to the opening of cracks in dilatation zones and to fluids in compression zones (Hillers et al., 2015; Renaud et al., 2014). Given these complexities, we therefore opt for using an average instead of conducting a station-by-station analysis. In addition, as using the peak value can distort the results due to edge effects generated nearby the fault. In fact, strain changes near the modeled fault's boundaries exhibit high strain anomalies that result from the numerical approximations used within the modeling. Therefore, to mitigate edge effects we opt for a median value of the highest 10% of dilatational strain changes, as suggested by Taira et al. (2018), yielding a 8.48 · 10⁻⁸ volumetric strain (Figure S8 in Supporting Information S1). This result implies that the magnitude of the localized velocity change (0.020%) is significantly larger than the strain. This discrepancy corroborates the hypothesis of non-linear elastic response of the shallow crust to low strain deformation. Indeed, for a linear elastic model, velocity changes are expected to be on the same order as the volumetric strain (Rivet et al., 2011).

A more quantitative evaluation of the relationship between strain and *dv/v* can be provided by the strain sensitivity parameter (Taira et al., 2018):
$$ \eta = \frac{dv/v}{\epsilon}, \qquad (9) $$

From Equation 9 we obtain a stress sensitivity η of roughly 2100, which is significantly larger than values obtained during coseismic shaking induced by large earthquakes. For example, the 2004 Mw 6 Parkfield earthquake (Chen et al., 2010) and the 2008 Mw 7.3 Wenchuan earthquake (Chen et al., 2010) have strain sensitivity values of 120 and 200, respectively.

Our estimated sensitivity also exceeds the values (from ~20 to ~60) obtained during seismic swarms observed in geothermal areas (Taira et al., 2018). On the other hand, comparable values of strain sensitivity (1000–1200) are observed in volcanic settings and are related to magma and fluid intrusion (Hirose et al., 2017; Takano et al., 2017).

The strain sensitivity is controlled by damage amount, or crack density and is known to rapidly decrease as function of depth, due to increment of confining pressure (Taira et al., 2018). Deeper layers of the crust are subjected to a higher confining pressure, which tends to counterbalance the opening and propagation of cracks. It is thus unlikely to observe a large strain sensitivity at seismogenic depth (3–5 km) as the one we report here (Brenguier, Shapiro, et al., 2008; Chen et al., 2010; Taira et al., 2018). Thus, the observed high value of strain sensitivity indicates a possible reduction in the confining pressure at depth, affecting the non-linear response of rocks (Johnson & Jia, 2005; Riviere et al., 2014). The presence of pressurized fluids at seismogenic depths can reduce confining pressure, promote crack formation and propagation and yield in return an increase of the strain sensitivity (Taira et al., 2018). We propose that the strong sensitivity observed is related to fractured and fluid permeable rocks in the shallow part of the ATF system.

This hypothesis agrees with observation of Piana Agostinetti et al. (2017), who proposed the occurrence of fluid-driven cracking of evaporitic rocks during the 2014 swarm. The complexity of the shallow fault system is also in agreement with the imbricated spatiotemporal evolution of seismicity observed by Essing and Poli (2022), during the 2014 swarm.

## 5. Conclusion

Our observation of significant velocity variation (0.033%) in response to small strain perturbation (~10⁻⁸) highlights the presence of complex rheology in the crust at seismogenic depth, which respond nonlinearly under small strain, and favor the occurrence of complex seismic swarms, rather than more classic mainshock aftershock earthquake sequences. Furthermore, we highlight how the monitoring of velocity variations, when integrated with seismological and geodetic observables, helps to track effective strain variations in faults, and provide quantitative information about the fault state. The same approach should be tested in other crustal environments affected by seismic swarms, particularly those where fluids are supposed to be involved.

## Data Availability Statement

The rain and the temperature data were obtained from the NASA Langley Research Center (LaRC) POWER Project funded through the NASA Earth Science/Applied Science Program (https://power.larc.nasa.gov/) (Stackhouse, 2021) and from a Joint Research Centre Data Catalog—EMO (Thiemig et al., 2022). The seismic data are available at Instituto Nazionale di Geofisica e Vulcanologia (http://iside.rm.ingv.it/instruments) (INGV, 2005). The GPS data (24 hr IGS14 final solutions) were downloaded from the Nevada Geodetic Laboratory website (http://geodesy.unr.edu) (Blewitt et al., 2018). The Coulomb software (Toda et al., 2011) was used to compute the volumetric strain associated with the 2013–2014 deformation episode.

## References

Agnew, D. C. (1981). Nonlinearity in rock: Evidence from Earth tides. *Journal of Geophysical Research*, 86(B5), 3969–3978. https://doi.org/10.1029/jb086ib05p03969

Bachura, M., Fischer, T., Doubravová, J., & Horálek, J. (2021). From earthquake swarm to a main shock–aftershocks: The 2018 activity in west bohemia/vogtland. *Geophysical Journal International*, 224(3), 1835–1848. https://doi.org/10.1093/gji/ggaa523

Barajas, A., Margerin, L., & Campillo, M. (2022). Coupled body and surface wave sensitivity kernels for coda-wave interferometry in a three-dimensional scalar scattering medium. *Geophysical Journal International*, 230(2), 1013–1029. https://doi.org/10.1093/gji/ggac091

## Acknowledgments

This research received funding from the European Research Council (ERC) under the European Union Horizon 2020 Research and Innovation Program (grant agreements 802777- MONIFAULTS). The authors thank GRICAD Grenoble Alpe Research - Scientific Computing and Data Infrastructure for their kind support during the computational part of this work. We also thank Satoshi Ide, Nori Nakata and Tim Clements for their constructive comments that significantly improved this manuscript.
Barajas, A., Poli, P., d'Agostino, N., Margerin, L., & Campillo, M. (2021). Separation of poroelastic and elastic processes of an aquifer from tectonic phenomena using geodetic, seismic, and meteorological data in the Pollino region, Italy. *Geochemistry, Geophysics, Geosystems*, 22(11), e2021GC009742. https://doi.org/10.1029/2021gc009742

Blewitt, G., Hammond, W., & Kreemer, C. (2018). Harnessing the GPS data explosion for interdisciplinary science [Dataset]. Eos, 99. https://doi.org/10.1029/2018eo104623

Boschelli, J., Moschetti, M. P., & Sens-Schönfelder, C. (2021). Temporal seismic velocity variations: Recovery following from the 2019 mw 7.1 Ridgecrest, California earthquake. *Journal of Geophysical Research: Solid Earth*, 126(4), e2020JB021465. https://doi.org/10.1029/2020jb021465

Bourouis, S., & Bernard, P. (2007). Evidence for coupled seismic and aseismic fault slip during water injection in the geothermal site of Soultz (France), and implications for seismogenic transients. *Geophysical Journal International*, 169(2), 723–732. https://doi.org/10.1111/j.1365-246x.2006.03325.x

Brenguier, F., Campillo, M., Hadzioioannou, C., Shapiro, N. M., Nadeau, R. M., & Larose, E. (2008). Postseismic relaxation along the San Andreas fault at Parkfield from continuous seismological observations. *Science*, 321(5895), 1478–1481. https://doi.org/10.1126/science.1160943

Brenguier, F., Campillo, M., Takeda, T., Aoki, Y., Shapiro, N., Briand, X., et al. (2014). Mapping pressurized volcanic fluids from induced crustal seismic velocity drops. *Science*, 345(6192), 80–82. https://doi.org/10.1126/science.1254073

Brenguier, F., Shapiro, N. M., Campillo, M., Ferrazzini, V., Duputel, Z., Coutant, O., & Nercessian, A. (2008). Towards forecasting volcanic eruptions using seismic noise. *Nature Geoscience*, 1(2), 126–130. https://doi.org/10.1038/ngeo104

Campillo, M., & Paul, A. (2003). Long-range correlations in the diffuse seismic coda. *Science*, 299(5606), 547–549. https://doi.org/10.1126/science.1078551

Cervelli, P., Murray, M., Segall, P., Aoki, Y., & Kato, T. (2001). Estimating source parameters from deformation data, with an application to the March 1997 earthquake swarm off the Izu peninsula, Japan. *Journal of Geophysical Research*, 106(B6), 11217–11237. https://doi.org/10.1029/2000jb900399

Chen, J. H., Froment, B., Liu, Q. Y., & Campillo, M. (2010). Distribution of seismic wave speed changes associated with the 12 May 2008 Mw 7.9 Wenchuan earthquake. *Geophysical Research Letters*, 37(18), L18302. https://doi.org/10.1029/2010gl044582

Daly, C., Neilson, R. P., & Phillips, D. L. (1994). A statistical-topographic model for mapping climatological precipitation over mountainous terrain. *Journal of Applied Meteorology and Climatology*, 33(2), 140–158. https://doi.org/10.1175/1520-0450(1994)033<0140:astmfm>2.0.co;2

De Barros, L., Cappa, F., Deschamps, A., & Dublanchet, P. (2020). Imbricated aseismic slip and fluid diffusion drive a seismic swarm in the Corinth Gulf, Greece. *Geophysical Research Letters*, 47(9), e2020GL087142. https://doi.org/10.1029/2020gl087142

De Barros, L., Guglielmi, Y., Rivet, D., Cappa, F., & Duboeuf, L. (2018). Seismicity and fault aseismic deformation caused by fluid injection in decametric in-situ experiments. *Comptes Rendus Geosciences*, 350(8), 464–475. https://doi.org/10.1016/j.crte.2018.08.002

De Gori, P., Lucente, F. P., & Chiarabba, C. (2015). Stressing of fault patch during seismic swarms in central Apennines, Italy. *Geophysical Research Letters*, 42(7), 2157–2163. https://doi.org/10.1002/2015gl063297

Essing, D., & Poli, P. (2022). Spatiotemporal evolution of the seismicity in the alto tiberina fault system revealed by a high-resolution template matching catalog. *Journal of Geophysical Research: Solid Earth*, 127(10), e2022JB024845. https://doi.org/10.1029/2022jb024845

Fiorillo, F. (2011). Tank-reservoir drainage as a simulation of the recession limb of karst spring hydrographs. *Hydrogeology Journal*, 19(5), 1009–1019. https://doi.org/10.1007/s10040-011-0737-y

Fischer, T., & Hainzl, S. (2021). The growth of earthquake clusters. *Frontiers in Earth Science*, 9, 638336. https://doi.org/10.3389/feart.2021.638336

Gualandi, A., Nichele, C., Serpelloni, E., Chiaraluce, L., Anderlini, L., Latorre, D., et al. (2017). Aseismic deformation associated with an earthquake swarm in the northern Apennines (Italy). *Geophysical Research Letters*, 44(15), 7706–7714. https://doi.org/10.1002/2017gl073687

Guglielmi, Y., Cappa, F., Avouac, J.-P., Henry, P., & Elsworth, D. (2015). Seismicity triggered by fluid injection-induced aseismic slip. *Science*, 348(6240), 1224–1226. https://doi.org/10.1126/science.aab0476

Hadzioioannou, C., Larose, E., Coutant, O., Roux, P., & Campillo, M. (2009). Stability of monitoring weak changes in multiply scattering media with ambient noise correlation: Laboratory experiments. *Journal of the Acoustical Society of America*, 125(6), 3688–3695. https://doi.org/10.1121/1.3125345

Hill, D. P. (1977). A model for earthquake swarms. *Journal of Geophysical Research*, 82(8), 1347–1352. https://doi.org/10.1029/jb082i008p01347

Hillers, G., Retailleau, L., Campillo, M., Inbal, A., Ampuero, J.-P., & Nishimura, T. (2015). In situ observations of velocity changes in response to tidal deformation from analysis of the high-frequency ambient wavefield. *Journal of Geophysical Research: Solid Earth*, 120(1), 210–225. https://doi.org/10.1002/2014jb011318

Hirose, T., Nakahara, H., & Nishimura, T. (2017). Combined use of repeated active shots and ambient noise to detect temporal changes in seismic velocity: Application to Sakurajima volcano, Japan. *Earth Planets and Space*, 69, 1–12. https://doi.org/10.1186/s40623-017-0613-7

Hobiger, M., Wegler, U., Shiomi, K., & Nakahara, H. (2012). Coseismic and postseismic elastic wave velocity variations caused by the 2008 Iwate-Miyagi Nairiku earthquake, Japan. *Journal of Geophysical Research*, 117(B9), B09313. https://doi.org/10.1029/2012jb009402

Hobiger, M., Wegler, U., Shiomi, K., & Nakahara, H. (2014). Single-station cross-correlation analysis of ambient seismic noise: Application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake. *Geophysical Journal International*, 198(1), 90–109. https://doi.org/10.1093/gji/ggu115

Hobiger, M., Wegler, U., Shiomi, K., & Nakahara, H. (2016). Coseismic and post-seismic velocity changes detected by passive image interferometry: Comparison of one great and five strong earthquakes in Japan. *Geophysical Journal International*, 205(2), 1053–1073. https://doi.org/10.1093/gji/ggw066

INGV. (2005). Rete sismica nazionale (RSN) [Dataset]. INGV. https://doi.org/10.13127/SD/X0FXNH7QFY

Johnson, P. A., & Jia, X. (2005). Nonlinear dynamics, granular media and dynamic earthquake triggering. *Nature*, 437(7060), 871–874. https://doi.org/10.1038/nature04015

Kato, A., Obara, K., Igarashi, T., Tsuruoka, H., Nakagawa, S., & Hirata, N. (2012). Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake. *Science*, 335(6069), 705–708. https://doi.org/10.1126/science.1215141

Li, G., & Ben-Zion, Y. (2023). Daily and seasonal variations of shallow seismic velocities in southern California from joint analysis of h/v ratios and autocorrelations of seismic waveforms. *Journal of Geophysical Research: Solid Earth*, 128(2), e2022JB025682. https://doi.org/10.1029/2022jb025682

Lobkis, O. I., & Weaver, R. L. (2003). Coda-wave interferometry in finite solids: Recovery of p-to-s conversion rates in an elastodynamic billiard. *Physical Review Letters*, 90(25), 254302. https://doi.org/10.1103/physrevlett.90.254302
Lohman, R., & McGuire, J. (2007). Earthquake swarms driven by aseismic creep in the Salton Trough, California. *Journal of Geophysical Research*, **112**(B4), B04405. https://doi.org/10.1029/2006jb004596

Maeda, T., Obara, K., & Yukutake, Y. (2010). Seismic velocity decrease and recovery related to earthquake swarms in a geothermal area. *Earth Planets and Space*, **62**(9), 685–691. https://doi.org/10.5047/eps.2010.08.006

Obermann, A., Froment, B., Campillo, M., Larose, E., Planès, T., Valette, B., et al. (2014). Seismic noise correlations to image structural and mechanical changes associated with the mw 7.9 2008 Wenchuan earthquake. *Journal of Geophysical Research: Solid Earth*, **119**(4), 3155–3168. https://doi.org/10.1002/2013jb010932

Obermann, A., Planès, T., Larose, E., Sens-Schönfelder, C., & Campillo, M. (2013). Depth sensitivity of seismic coda waves to velocity perturbations in an elastic heterogeneous medium. *Geophysical Journal International*, **194**(1), 372–382. https://doi.org/10.1093/gji/ggt043

Okada, Y. (1992). Internal deformation due to shear and tensile faults in a half-space. *Bulletin of the Seismological Society of America*, **82**(2), 1018–1040. https://doi.org/10.1785/bssa0820021018

Ostrovsky, L., & Johnson, P. A. (2001). Dynamic nonlinear elasticity in geomaterials. *La Rivista del Nuovo Cimento*, **24**(7), 1–46. https://doi.org/10.1007/bf03548898

Passarelli, L., Rivalta, E., Jónsson, S., Hensch, M., Metzger, S., Jakobsdóttir, S. S., et al. (2018). Scaling and spatial complementarity of tectonic earthquake swarms. *Earth and Planetary Science Letters*, **482**, 62–70. https://doi.org/10.1016/j.epsl.2017.10.052

Piana Agostinetti, N., Giacomuzzi, G., & Chiarabba, C. (2017). Seismic swarms and diffuse fracturing within Triassic evaporites fed by deep degassing along the low-angle alto tiberina normal fault (central Apennines, Italy). *Journal of Geophysical Research: Solid Earth*, **122**(1), 308–331. https://doi.org/10.1002/2016jb013295

Poli, P., Cabrera, L., Flores, M., Báez, J., Ammirati, J., Vásquez, J., & Ruiz, S. (2022). Volcanic origin of a long-lived swarm in the central bransfield basin, Antarctica. *Geophysical Research Letters*, **49**(1), e2021GL095447. https://doi.org/10.1029/2021gl095447

Poli, P., Marguin, V., Wang, Q., d'Agostino, N., & Johnson, P. (2020). Seasonal and coseismic velocity variation in the region of L'Aquila from single station measurements and implications for crustal rheology. *Journal of Geophysical Research: Solid Earth*, **125**(7), e2019JB019316. https://doi.org/10.1029/2019jb019316

Poupinet, G., Ellsworth, W., & Frechet, J. (1984). Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California. *Journal of Geophysical Research*, **89**(B7), 5719–5731. https://doi.org/10.1029/jb089ib07p05719

Renaud, G., Rivière, J., Larmat, C., Rutledge, J., Lee, R., Guyer, R., et al. (2014). In situ characterization of shallow elastic nonlinear parameters with dynamic acoustoelastic testing. *Journal of Geophysical Research: Solid Earth*, **119**(9), 6907–6923. https://doi.org/10.1002/2013jb010625

Richter, T., Sens-Schönfelder, C., Kind, R., & Asch, G. (2014). Comprehensive observation and modeling of earthquake and temperature-related seismic velocity changes in northern Chile with passive image interferometry. *Journal of Geophysical Research: Solid Earth*, **119**(6), 4747–4765. https://doi.org/10.1002/2013jb010695

Rivet, D., Brenguier, F., & Cappa, F. (2015). Improved detection of preruptive seismic velocity drops at the Piton de La Fournaise volcano. *Geophysical Research Letters*, **42**(15), 6332–6339. https://doi.org/10.1002/2015gl064835

Rivet, D., Campillo, M., Shapiro, N. M., Cruz-Atienza, V., Radiguet, M., Cotte, N., & Kostoglodov, V. (2011). Seismic evidence of nonlinear crustal deformation during a large slow slip event in Mexico. *Geophysical Research Letters*, **38**(8), L08308. https://doi.org/10.1029/2011gl047151

Riviere, J., Remillieux, M. C., Ohara, Y., Anderson, B., Haupert, S., Ulrich, T., & Johnson, P. A. (2014). Dynamic acousto-elasticity in a fatigue-cracked sample. *Journal of Nondestructive Evaluation*, **33**(2), 216–225. https://doi.org/10.1007/s10921-014-0225-0

Roland, E., & McGuire, J. J. (2009). Earthquake swarms on transform faults. *Geophysical Journal International*, **178**(3), 1677–1690. https://doi.org/10.1111/j.1365-246x.2009.04214.x

Scuderì, M., Marone, C., Tinti, E., Di Stefano, G., & Collettini, C. (2016). Precursory changes in seismic velocity for the spectrum of earthquake failure modes. *Nature Geoscience*, **9**(9), 695–700. https://doi.org/10.1038/ngeo2775

Sens-Schönfelder, C., & Wegler, U. (2006). Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia. *Geophysical Research Letters*, **33**(21), L21302. https://doi.org/10.1029/2006gl027797

Sens-Schönfelder, C., & Wegler, U. (2011). Passive image interferometry for monitoring crustal changes with ambient seismic noise. *Comptes Rendus Geoscience*, **343**(8–9), 639–651. https://doi.org/10.1016/j.crte.2011.02.005

Shapiro, N. M., & Campillo, M. (2004). Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise. *Geophysical Research Letters*, **31**(7), L07614. https://doi.org/10.1029/2004gl019491

Shapiro, N. M., Campillo, M., Stehly, L., & Ritzwoller, M. H. (2005). High-resolution surface-wave tomography from ambient seismic noise. *Science*, **307**(5715), 1615–1618. https://doi.org/10.1126/science.1108339

Shapiro, S. A., Huenges, E., & Borm, G. (1997). Estimating the crust permeability from fluid-injection-induced seismic emission at the KTB site. *Geophysical Journal International*, **131**(2), F15–F18. https://doi.org/10.1111/j.1365-246x.1997.tb01215.x

Sheng, Y., Mordret, A., Sager, K., Brenguier, F., Boué, P., Rousset, B., et al. (2022). Monitoring seismic velocity changes across the San Jacinto fault using train-generated seismic tremors. *Geophysical Research Letters*, **49**(19), e2022GL098509. https://doi.org/10.1029/2022gl098509

Snieder, R., Grét, A., Douma, H., & Scales, J. (2002). Coda wave interferometry for estimating nonlinear behavior in seismic velocity. *Science*, **295**(5563), 2253–2255. https://doi.org/10.1126/science.1070015

Stackhouse, P. (2021). NASA Langley Research Center (LARC) power project [Dataset]. NASA. Retrieved from https://power.larc.nasa.gov/

Taira, T., Nayak, A., Brenguier, F., & Manga, M. (2018). Monitoring reservoir response to earthquakes and fluid extraction, Salton Sea geothermal field, California. *Science Advances*, **4**(1), e1701536. https://doi.org/10.1126/sciadv.1701536

Takano, T., Nishimura, T., & Nakahara, H. (2017). Seismic velocity changes concentrated at the shallow structure as inferred from correlation analyses of ambient noise during volcano deformation at Izu-Oshima, Japan. *Journal of Geophysical Research: Solid Earth*, **122**(8), 6721–6736. https://doi.org/10.1002/2017jb014340

Thiemig, V., Gomes, G. N., Skøien, J. O., Ziese, M., Rauthe-Schöch, A., Rustemeier, E., et al. (2022). EMO-5: A high-resolution multi-variable gridded meteorological dataset for Europe [Dataset]. Earth System Science Data, **14**, 3249–3272. https://doi.org/10.5194/essd-14-3249-2022

Toda, S., Stein, R. S., Sevilgen, V., & Lin, J. (2011). Coulomb 3.3 graphic-rich deformation and stress-change software for earthquake, tectonic, and volcano research and teaching—User guide [Software]. *US Geological Survey Open-File Report*, **1060**, 63. Retrieved from https://pubs.usgs.gov/of/2011/1060/

Vadacca, L., Casarotti, E., Chiaraluce, L., & Cocco, M. (2016). On the mechanical behaviour of a low-angle normal fault: The Alto Tiberina fault (northern Apennines, Italy) system case study. *Solid Earth*, **7**(6), 1537–1549. https://doi.org/10.5194/se-7-1537-2016

Vidale, J. E., & Shearer, P. M. (2006). A survey of 71 earthquake bursts across southern California: Exploring the role of pore fluid pressure fluctuations and aseismic slip as drivers. *Journal of Geophysical Research*, **111**(B5), B05312. https://doi.org/10.1029/2005jb004034

Wang, Q.-Y., Brenguier, F., Campillo, M., Lecointre, A., Takeda, T., & Aoki, Y. (2017). Seasonal crustal seismic velocity changes throughout Japan. *Journal of Geophysical Research: Solid Earth*, **122**(10), 7987–8002. https://doi.org/10.1002/2017jb014307
Wegler, U., & Sens-Schönfelder, C. (2007). Fault zone monitoring with passive image interferometry. *Geophysical Journal International*, *168*(3), 1029–1033. https://doi.org/10.1111/j.1365-246x.2006.03284.x

Yamamura, K., Sano, O., Utada, H., Takei, Y., Nakao, S., & Fukao, Y. (2003). Long-term observation of in situ seismic velocity and attenuation. *Journal of Geophysical Research*, *108*(B6), 2317. https://doi.org/10.1029/2002jb002005

Zhao, P., Peng, Z., & Sabra, K. G. (2010). Detecting remotely triggered temporal changes around the Parkfield section of the San Andreas fault. *Earthquake Science*, *23*(5), 497–509. https://doi.org/10.1007/s11589-010-0748-0