The Coda Normalization Method, introduced by Aki (1980), is a technique used to calculate the attenuation of the S-wave from waveforms recorded at multiple seismic stations. This method has been widely applied in both volcanic (De Siena et al., 2014; Castro-Melgar et al., 2021; Bala et al., 2024 and references therein) and tectonic zones (Sedaghti and Pezeshk, 2016; Sketsiou et al., 2021; Toker and Şahin, 2022; Talone et al., 2023, and references therein). Moreover, defining the attenuation characteristics of an area turns out to be fundamental in the seismic hazard assessment (Lombardi et al., 2005; Akinci et al., 2010). This method consists of the ratio between the direct S-wave amplitude and the amplitude of the coda waves, effectively removing the site, source intensity and instrumental effects. By using this technique, we can also retrieve the total quality factor and its frequency dependence, Q(f), in the area. The attenuation of the S wave can then be calculated as ln A S f , r r − γ R ϕ θ A c f , t c = s − f 2 Q β V s r (3) where A_s(f, r) is the S wave spectral amplitude, A_c(f, t_c) is the coda wave spectral amplitude, f is the frequency, r is the source-receiver distance, R_θφ is the source radiation pattern (θ and φ the azimuth and take-off angle for a source–receive ray), γ the geometrical spreading exponent, t_c is the coda wave lapse time, s is a function of the source of S and coda waves, and site effect at the receiver, Q_β is the S-wave quality factor and V_S the average S-wave velocity in the medium (here, 3.5 km/s). The slope of the least-squares fit of this equation to the corresponding observational data then determines the total attenuation Q_β. The starting lapse time for the coda window is set at 75 s, which is more than twice the S wave arrival, while the body window length is set at 10 s (Figure 2). This extended lapse time and body window are chosen to accommodate the large hypocentral distances (up to 100 km) in our dataset, ensuring that the energy peak is captured within the body window and preventing contamination of the coda window with the S wave energy. For source-receiver distances up to 100 km, the S-wave phase arrives around 28 s from the origin, with the coda starting at 42 s, providing over 30 s of signal time for the transition from single scattering to a diffusion regime.

We applied a fourth-order band-pass Butterworth filter, forward and backward, in four frequency bands (1–2 Hz, 2–4 Hz, 4–8 Hz, and 8–16 Hz), centered at frequencies of 1.5, 3, 6, and 12 Hz (Gabrielli et al., 2022). Figure 3 presents an example of the pre-sequence phase coda normalization analysis across these four frequency bands. To compute the envelopes, we used a Hilbert transform and then smoothed them with a moving window of duration equal to the inverse of eight times the central frequency. This widely used procedure, applied in many studies on similar topics, facilitates useful comparisons across comparable results worldwide (Calvet et al., 2013; Borleanu et al., 2017; Sketsiou et al., 2021).

An initial analysis of the entire dataset, including near-source waveforms, highlighted a notable effect within a hypocentral distance of 0–20 km, where a significant increase in the energy ratio caused a bias in the calculation of Q_β at lower frequencies (1.5 and 3 Hz) (Figure 4). This behavior at close distances (0–20 km) may be influenced by factors such as scattering due to complex heterogeneities and an azimuthal component of geometrical spreading. Given the difficulty in discerning the dominant effect within this 0–20 km range, we chose to exclude waveforms within the first 20 km. Consequently, for the Q_β estimation, we considered only seismic waves with hypocentral distance between 20 and 100 km. We also adjusted the geometrical spreading values, as displayed in Table 3. This adjustment helped reduce the impact of unusually high energy values at a low (0–30 km) hypocentral distance.

Because the coda normalization method provides a single value of total attenuation for each dataset, it does not facilitate 2D mapping. To visualize spatial variations, we used the coda quality factor (Q_c) as a proxy for intrinsic attenuation.

Coda waves, which appear as the tail end of a seismogram after the S-wave arrival, result from the complexities seismic waves encounter while traveling through the Earth’s crust (Aki and Chouet, 1975). A key characteristic of coda waves is their attenuation, where the wave envelope diminishes over time from the event’s origin (Zhang et al., 2021).

Various factors influence coda wave behavior, including single scattering, multiple scattering, and diffusion, depending on frequency and time-lapse (Sato et al., 2012). Shapiro et al. (2000) focused on the diffusion approximation of multiple scattering models and demonstrated that late coda radiation tends to be a diffusive process. In their work, they analyzed how coda waves transition into a diffusive regime. In this context, diffusion represents the complete redistribution of wave energy throughout the medium. In the case of small intrinsic dissipation (high intrinsic-Q), the radiation emitted by the source is not affected by a substantial loss of energy.

Calvet and Margerin (2013) expanded on the understanding of coda waves by studying the energy decay and scattering in seismology, focusing on both diffusion and absorption phenomena. They highlighted how coda waves can be dominated by single scattering in the early coda, and then exhibit a diffusive behavior in the late coda. They also made comparisons between scattering-induced energy losses and intrinsic absorption, emphasizing the different physical processes governing each phenomenon.

In extreme synthesis, for late coda, the coda energy decay (Q-coda coefficient, Q_c, as explained in the following) tends to the value of intrinsic Q, Q_i.

The decay rate of the coda envelope is quantified by the inverse coda quality factor Q⁻¹ _c (Aki and Chouet, 1975), and it is expressed by the equation: E t , f = S f t − α ⁡ exp − 2 π f t Q c (4)

Where E (t, f) represents the power spectral density, S(f) the frequency-dependent source-site term, t the time-lapse from the origin, and α is equal to 3/2 when propagation occurs in a single layer defined by a multiple scattering regime (Paasschens, 1997; Calvet et al., 2013).

Q⁻¹ _c is derived by linearizing Equation 4 through straight-line fitting, as in Equation 5: ln E t , f ∙ t α 2 π f = ln S f 2 π f − 1 Q c t (5)

For Q⁻¹ _c analysis, the coda window starts at 2 times the S wave arrival and has a length of 20 s, following the tests done by Gabrielli et al. (2022) (Figure 2).

After determining the inverse quality factor Q⁻¹ _c at each central frequency (1.5–3–6–12 Hz) for all station-source pairs, we mapped the spatial variations in Q⁻¹ _c using kernel functions applicable to the diffusive regime (Del Pezzo et al., 2016). These kernels, related to the diffusive processes, provide higher spatial sensitivity in heterogeneous media, enhancing the illumination around the seismic path. This approach has been applied in both volcanic (De Siena et al., 2017; Gabrielli et al., 2020) and tectonic (Napolitano et al., 2020; Borleanu et al., 2023) settings.

After applying the kernels, we mapped the Q⁻¹ _c using a regionalization approach, which involves calculating a weighted average of Q_c from all the rays passing through each block. Although this regionalization method cannot include a resolution test, like a checkerboard pattern, due to its non-inversion nature, we evaluated the validity of the results using a hit count map (Calvet et al., 2013; Gabrielli et al., 2022). This approach, widely used in attenuation studies, serves as a proxy for resolution testing (Takahashi et al., 2007; Calvet et al., 2013; Gabrielli et al., 2022; Gabrielli et al., 2023; Napolitano et al., 2023). In Section 5.2, we are going to present the hit count map alongside the results, incorporating only cells crossed by more than 10 rays, in accordance with the approach outlined by Calvet et al. (2013).

To present the robustness and reliability of the spatial variation of attenuation, we employed the technique developed by Del Pezzo and Ibáñez (2020) to show the spatial distribution of the normalized standard deviation σ associated with the single Q values. This new approach involves defining “Resolution” as the inverse of the normalized standard deviation of the weighted average computed for each grid cell. The weights are determined by the numerical value of the kernel function at the center of each cell. High values of resolution are associated with the more reliable estimates in each cell, and vice versa. Del Pezzo and Ibáñez (2020) and Castro-Melgar et al. (2021) applied this new representation of uncertainty in the volcanic area of the Aeolian Islands. This kind of test can be considered the most suitable for the projection mapping of the weighting function of Q_c that we performed. Here, for the first time, we applied this technique in tectonic regions with kernels having longer distances than the one used for volcanic zones (>20 km). We furthermore tested the applicability of these kernels, as presented in Supplementary Figure S1 provided in Supplementary Material.

We performed a study of Coda Normalization for the two main sequences (the pre-sequence and the 2016-2017 sequence) as well as for the three sub-sequences of Amatrice, Visso, and Norcia (Figure 5). In Figure 5, we plotted the distribution of the Energy ratios derived from Equation 3 against the hypocentral distance for each sequence, and at the central frequency of 1.5 Hz. Considering that the coda normalization has been performed also for the frequencies of 3, 6, and 12 Hz, we are able to obtain the quality factor Q and its frequency dependence, following Q f = Q 0 f n , where Q₀ is the quality factor Q at 1 Hz. It is important to mention that, although in Figure 5, the Visso phase appears smaller than the other sequences, the trend as a function of distance is still clearly observable. So, the consistency of the observed trend (which is marked by a decrease in the energies ratio with increasing distances), even with fewer data points, supports the conclusion of a lower value of Q (with respect to the values calculated for the other sequences), independently from the dataset size.

The frequency dependence of attenuation for the pre-sequence is described by Q_β(f) = 93 ± 25 f ^0.58±0.11, while during the sequence, it is Q_β(f) = 80 ± 22 f ^0.66±0.15. For specific sequences of Amatrice, Visso, and Norcia, the frequency dependences are: Amatrice: Q_β(f) =98 ± 25 f ^0.56±0.14; Visso: Q_β(f) =58 ± 21f ^0.8±0.16; Norcia: Q_β(f) =74 ± 20f ^0.7±0.15. These results are summarized in Table 4.

These results suggest that the crustal S-wave attenuation and its frequency dependence varied significantly during the main events of the AVN sequence. Between the pre-sequence and the sequence, Q₀ values decreased by 14%, indicating an increase in attenuation during this period. Throughout the sequence, particularly during the Visso phase (where Q₀ = 58, from Q_β(f) = Q₀fⁿ), there is an apparent decrease in Q_β, with a 27% decrease compared to the sequence and a 37% decrease compared to the pre-sequence. This decrease may be associated with changes in fluid distribution during the sequence, as observed in other attenuation studies (Gabrielli et al., 2022; Castro et al., 2022; Akinci et al., 2022). Gabrielli et al. (2022) observed, in their Figure 9, possible fluid migration during the Visso phase by analyzing the plot of the average Q_c ⁻¹ at two different stations. At the northernmost station located in the Visso area, we detected a 40% increase in attenuation just after the pre-sequence and during the Visso and Norcia phases, comparable to the 37% variation we found with the coda normalization. It is important to note that the attenuation represented by Gabrielli et al. (2022) is Q_c, not the Q_β, even if both can be influenced by the fluid content. This fluid migration towards the Visso area was also identified by Malagnini et al. (2022), who calculated fluid flow patterns and suggested a northward fluid movement after the Amatrice mainshock.

Alongside Q_β, there is a corresponding variation in the frequency parameter, n, which increases as Q₀ decreases and vice versa. The observed variations in Q_β values are comparable with those reported by Akinci et al. (2022) using the random vibration theory (RVT) technique (Tables 1, 4). In comparison to other studies (Table 1), the Q values from the GIT (Generalized Inversion Technique) analysis by Morasca et al. (2023 - Q = 247) and Pacor et al. (2016 - Q = 290) are higher than those obtained using the coda normalization used in this work (see Tables 1, 4). However, their n values are way lower than those reported in this study.

In this section, we present the spatial variation of coda waves’ attenuation in the study area. Figures 6A, 7 represent areas of high attenuation Q_c ⁻¹ in dark orange and areas of low attenuation Q_c ⁻¹ in blue. The results are plotted after selecting cells based on the ray count per cell, considering only those crossed by more than 10 rays (hit count map, Figure 6B). The interpretation is further supported by the standard deviation σ map, which presents the relative cell resolution on calculated attenuation values (Figure 6). As described in the Methodology section, we combined the pre-sequence and sequence datasets to examine coda wave attenuation variation across the entire Central Apennines area rather than focusing solely on the 2016–2017 seismic sequence (Figure 6A). Additionally, we provide maps for the individual sequences of Amatrice, Visso, and Norcia (Figure 7). We focus on 1.5 Hz as it allows us to observe the primary structures and dynamics of the tectonic setting, enabling us to visualize structures with a linear dimension on the order of 2-3 km. The results at the higher frequency of 12 Hz are shown in Supplementary Figure S2 provided in Supplementary Material.

In the center of the map (Figure 6A), we observe that high attenuation anomalies are located in the same regions identified in previous studies on coda attenuation (Gabrielli et al., 2022). The expanded study area, the introduction of new data, and the regionalization approach have also revealed significant attenuation along the Adriatic coast, likely due to the presence of foredeep deposits in the area. Additionally, the high attenuation anomaly extends along the thrust front of the Monti Sibillini (MST, Figure 1), consistent with findings from other attenuation studies (Gabrielli et al., 2022; Colavitti et al., 2022; Morasca et al., 2023). We identify a similar trend at a higher frequency of 12 Hz (Supplementary Figure S2). The Monti Sibillini thrust is considered an important element in the origin of the AVN seismic sequence, as seen also by the fieldwork studies of Curzi et al. (2024), where they suggest that this thrust may represent a barrier for the deep fluid migration at the hypocentral depth of the present-day earthquake. The involvement of this thrust in the seismic sequence starting with the Amatrice mainshock, has been recognized also by the aftershock distribution during the seismic sequence, as seen by Bonini et al. (2019) and Del Ventisette et al. (2021). This high attenuation trend also extends northward and westward to the Perugia region, an area that has been subject to a seismic sequence in 1997–1998 (Chiaraluce et al., 2004).

Moreover, the high attenuation anomaly extending northward, partially visible during the Norcia phase (Figure 7), aligns with observations of crustal velocity variations during the 2016-2017 Amatrice-Visso-Norcia (AVN) sequence reported by Soldati et al. (2019). Using noise cross-correlation techniques, Soldati et al. identified a sudden drop in seismic velocity following each mainshock, with the most significant fluctuations occurring along the NE direction. These variations extended toward and along the Adriatic coast, particularly in sandstone-rich foredeep deposits. This phenomenon is likely driven by the porous and hydrated nature of these sediments, which respond to seismic wave propagation by increasing pore pressure. Elevated pore pressure reduces the effective stress within the sediments, causing a decrease in seismic velocity. Such processes underscore the dynamic role of fluid-saturated geological formations in modulating seismic wave behavior during an active sequence. Our findings are consistent with these results, as demonstrated by the observed increase in attenuation during the AVN sequence. Specifically, coda normalization analysis reveals a 24% decrease in Q (indicating an increase in attenuation) between the Amatrice and Norcia mainshocks (Table 4). This change points to a rise in fluid content in the area between the Visso and Norcia events, further supporting the hypothesis of fluid-driven processes. The Q_c⁻¹ mapping (Figures 6A, 7) also highlights spatially correlated attenuation anomalies, emphasizing the impact of fluid redistribution and sediment properties on seismic wave propagation.

In southern Central Italy, we observe a low attenuation zone extending east of Rome. This low Q_c ⁻¹ zone is likely associated with changes from high to medium CO₂ flux and high to low conductive heat flow, as noted by Di Luccio et al. (2022), Figures 3A, B. It is interesting to see that the decrease in attenuation aligns with the variations in CO₂ and conductive heat seen by Di Luccio et al. (2022), following an almost NE-SW direction. These variations are influenced by the Apennines structure, which separates the western Tyrrhenian side, characterized by diffusive mantle-source CO₂ release, degassing, and volcanic activity, from the Adriatic foredeep to the east, where no CO₂ emission is observed (Chiodini et al., 2020; Di Luccio et al., 2022). The lower attenuation area in Figure 6A aligns with the edge of these variations, marking the transition between volcanic gas emissions and the sedimentary domain of the Apennine chain.

The high attenuation anomaly in Figure 6A, observed in the foredeep of the Adriatic Coast (south of Ancona), corresponds with findings from Carannante et al. (2013), who identified a high Vp/Vs ratio in the same area between the depths of 0 and 4 km, suggesting the presence of highly hydrated Plio-Pleistocene sediments (Scrocca et al., 2007). At a greater depth of 12 km, Carannante et al. (2013) detected a contrast in Vp/Vs ratios, indicating the separation of two different flexing portions of the Adria lithosphere caused by a raised area of the basement. It is important to mention that while Carannante et al. (2013) conducted a 3D tomography, our study is based on 2D imaging, so we cannot determine which possible interpretation among those features can be the cause of our high attenuation anomaly.

To examine the attenuation changes from 2011 to 2017, we calculated the percentage difference (Q_c ⁻¹ variation) between the Q_c ⁻¹ of the 2016 sequence and the pre-sequence phase (Figure 8). Imaging of the individual Pre-Sequence and Sequence phases is provided in Supplementary Figure S3. It is important to note that the coverage between the two sequences is similar, but not equal, as the Adriatic coast is partially more covered during the pre-sequence phase, due to the new events added in the dataset.

The main variations are observed to the west of the AVN mainshocks and of the Monti Sibillini thrust, with changes in Q_c ⁻¹ of approximately 20–30%. These differences are comparable with the variation in coda-wave attenuation and scattering reported by Gabrielli et al. (2022); Gabrielli et al. (2023), though the additional seismicity in our dataset reveals further differences. Indeed, these variations are located where the seismic events of 2016 occurred, recording an increment of attenuation with respect to the pre-sequence. These variations appear as blue anomalies along the foredeep and Adriatic Coast, indicating visible negative variation (−40%), as well as changes in the northern part of the study area.

We also calculated an average of Q_c and its frequency dependence based on the spatial distribution (Table 5) ( Q c f = Q c 0 f n ). This comparison (Table 5; Figure 9) reveals time and space variations in Q_c0 and Q₀. While the frequency dependence of Q_c0 remains relatively constant across the four seismic phases, more substantial variations are observed in the frequency pattern of Q₀ (Figure 9). Q₀ provides an indirect estimate of total Q at 1 Hz, while Qc₀ is an estimate of Q_c at 1 Hz. The differences in Q₀ frequency pattern suggest potential changes in scattering attenuation (Q_sc) during the seismic phases, as noted by Gabrielli et al. (2023). Gabrielli et al. (2022) reported both temporal and spatial variations in Q_c, noting changes in 1/Q_c at two stations, one located north and the other south of the seismic sequence (see Figures 7, 9 in Gabrielli et al., 2022). These fluctuations are likely due to fluid migration within the system, as inferred also from the attenuation mapping of the sequences. However, in Gabrielli et al. (2022), the frequency-dependent 1/Q_c remained stable between the pre-sequence (133f ^0.75) and the sequence (128f ^0.75), which are close values to the one found in this study (102 ± 13f ^0.67±0.05 for the pre-sequence and 114 ± 17f ^0.69±0.06 for the sequence). It is important to note that although both papers share the same dataset, in Gabrielli et al. (2022), the dataset-particularly for the pre-sequence- was reduced for tomographic purposes. In contrast, in this study, through the application of regionalization, we utilized nearly the entire dataset and included additional events, which contributed to the observed variations in Q values and their frequency dependence.