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

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.


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 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 weakenning 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 asses rocks characteristics including the dynamic moduli such as shear, bulk, and Young's modulus, the mechanical hysteresis, the resonnance 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 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.

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 aut-ocorrelation (AC) and the single station crosscorrelation (SC), yield an estimation of the EGF (Hobiger et al., 2014(Hobiger et al., , 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 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.

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.

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

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

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

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-avergaged 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: 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: If we discretize Equation 3 with a first order approximation we end up with the following expression: Since our time resolution is of one day, we can reduce Equation 4 to: 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 4 .We later assess the model fit to the measured data with the misfit σ(k) 2 (6).The best fit corresponds to the k value that gives the minimal misfit.
Journal of Geophysical Research: Solid Earth 10.1029/2023JB028232 MIKHAEL ET AL.
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.

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): 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 simoultaneous 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 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 show 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 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 indentify 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).Additionaly, 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.

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 seismogenic 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 w 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: with d being the amplitude of the velocity drop associated with the seismic swarm and D r (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 4 to 1.38 ⋅ 10 4 .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 entrappement 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 whithin 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 8 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): From Equation 9we 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 M w 6 Parkfield earthquake (Chen et al., 2010) and the 2008 M w 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 fluiddriven 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.

Conclusion
Our observation of significant velocity variation (0.033%) in response to small strain perturbation (∼10 8 ) 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.

Figure 1 .
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.

Figure 2 .
Figure 2. Stack of relative velocity variations measured from autocorrelation 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.

Figure 3 .
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.

Figure 4 .
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.

Figure 5 .
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 corespond to time periods exhibiting large velocity drops.(b) Number of missing stations from the dv/v averaging process in gray dots.

Figure 6 .
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 correponding 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.

Table 1 σ
(k) 2 Values Between Synthetic and Observed dv/v for All Coda Lapse Time Windows and Each of the Three Models Used 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.
MIKHAEL ET AL.