Fast Changes in Seismic Attenuation of the Upper Crust due to Fracturing and Fluid Migration: The 2016–2017 Central Italy Seismic Sequence

The Amatrice–Visso–Norcia seismic sequence struck Central Italy across the Apenninic normal fault system in 2016. Fluids likely triggered the sequence and reduced the stability of the fault network following the first earthquake (Amatrice, Mw 6.0), with their migration nucleating the Visso (Mw 5.9) and Norcia (Mw 6.5) mainshocks. However, both spatial extent and mechanisms of fluid migration and diffusion through the network remain unclear. High fluid content, enhanced permeability, and pervasive microcracking increase seismic attenuation, but different processes contribute to different attenuation mechanisms. Here, we measured and mapped peak delay time and coda attenuation, using them as proxies of seismic scattering and absorption before and during the sequence. We observed that the structural discontinuities and lithology control the scattering losses at all frequencies, with the highest scattering delineating carbonate formations within the Gran Sasso massif. The Monti Sibillini thrust marks the strongest contrasts in scattering, indicating a barrier for northward fracture propagation. Absorption does not show any sensitivity to the presence of these main geological structures. Before the sequence, low-frequency high-absorption anomalies distribute around the NW-SE-oriented Apennine Mountain chain. During the sequence, a high-absorption anomaly develops from SSE to NNW across the seismogenic zone but remains bounded north by the Monti Sibillini thrust. We attribute this spatial expansion to the deep migration of CO2-bearing fluids across the strike of the fault network from a deep source of trapped CO2 close to the Amatrice earthquake. Fluids expand SSE-NNW primarily during the Visso sequence and then diffuse across the fault zones during the Norcia sequence.


INTRODUCTION
Between August and October 2016, 25.600 earthquakes (0.1 < M < 6.5) struck the Central Apennine chain in Central Italy (Chiaraluce et al., 2017). Eight events had M w > 5.0, with three mainshocks occurring near Amatrice (M w 6.0, August 24, 2016), Visso (M w 5.9, October 26, 2016), and Norcia (M w 6.5, October 30, 2016). The Amatrice-Visso-Norcia sequence (AVN) developed on the NNW-SSE-trending normal fault systems of Mt. Vettore (north) and Mt. Gorzano (south), activating a 70 km × 10 km area trending parallel to the major axis of the Central Apennines (Improta et al., 2019) (Figure 1). The fault mechanisms were consistent with the extension of the Central Apennines, which coexists with the opening of the Tyrrhenian back-arc basin (Chiaraluce et al., 2017). Digital seismic networks were deployed immediately after the Amatrice earthquake, recording events with M > 0.1 and making the AVN the best-monitored earthquake sequence striking Italy. The AVN has been modeled as a cascading rupture triggered by fluids entering the fault systems (Walters et al., 2018). Multiple mainshock sequences are common features in the Apennines (Improta et al., 2019), as the 1980 M s 6.9 of Irpinia (Bernard and Zollo, 1989), the 1997 M w 6.0 of Umbria-Marche , the 2009 M w 6.1 of L'Aquila (Valoroso et al., 2013), and the 2012 M L 5.0 of Pollino (De Matteis et al., 2021). This phenomenon is a consequence of the complex interactions between adjacent fault systems and triggering processes that include static and dynamic stress transfer between fault systems, pore pressure variations associated with fluid migration (Ventura and Di Giovambattista, 2013), and co-seismic release of a deep source of trapped CO 2 (Miller et al., 2004;Chiarabba et al., 2009a;Malagnini et al., 2012).
Waves attenuate due to scattering and absorption (Sato et al., 2012). The increase of seismic absorption relative to scattering losses and the frequency dependence of the attenuation measurements both support the inference of fluids permeating faults during the AVN (Akinci et al., 2020). The seismic waves yield energy when traveling through complex fluid-filled fault networks (Sketsiou et al., 2020). Seismic attenuation tomography has thus the potential to infer the effect of increased fluid content and migration. Chiarabba et al. (2009a) used P-and S-wave attenuation to track fluid migrations across the Central-Northern Apennine fault networks, identifying them as the driving mechanism of the cause of the 1997 Umbria-Marche seismic sequence (M w 5.7 and 6.0). Sketsiou et al. (2020) imaged total attenuation across the Pollino range (Southern Italy), detecting two fluid reservoirs hosting fluids that laterally migrated at depth, producing the 2010-2014 Pollino seismic sequence (M L 4.3 and 5.0- Brozzetti et al. (2017) and references therein).
Scattering and absorption can be mapped in space using two seismic attributes: FIGURE 1 | Simplified geological map of Central Italy (after Di Bucci et al., 2021). Blue: carbonates of the Umbria-Marche (UMD) and Lazio-Abruzzi (LAD) domains; yellow: foredeep domain (LF, Laga Formation). The seismic network used for the sequence period is represented with red triangles. Yellow stars correspond to the Amatrice (M w 6), Visso (M w 5.9), and Norcia (M w 6.5) mainshocks. The black and purple boxes are the Mt. Gorzano (Amatrice) and Mt. Vettore (Norcia) fault planes, respectively. Blue circles are the Nuria-Velino-Giano (NVG) and Gran Sasso (GS) aquifers. The study area within the Italian Peninsula framework is shown in the inset on the right panel.
1) The peak delay time, defined as the time difference between the S-wave arrival and the maximum amplitude of the event (Takahashi et al., 2007); 2) The attenuation of coda waves (Q c −1 ). Coda waves are wavetrains coming after the S-wave packet, and Q c quantifies the decay rate of the coda envelope with increasing lapse time from the origin time of the earthquake (Aki and Chouet, 1975).
At crustal-scale and in the far-field, the peak delay measures multiple forward scattering due to random inhomogeneities (Takahashi et al., 2007;. The Markov approximation (Saito et al., 2002) can model its variations in space and frequency at frequencies higher than 1.0 Hz. Takahashi et al. (2007) first used peak delay times to map volcanic areas in northeast Japan. High peak delays mark fractured volumes across the Pyrenees , the eastern portion of the Siletz terrane in the Western US (De Siena et al., 2016), and the Vrancea region in Romania (Borleanu et al., 2017). Within fault networks, peak delay increases are visible at all frequencies when waves cross fractured geological volumes, as the carbonates of the Pollino seismic gap, in Southern Italy .
Single scattering, multiple scattering, and diffusion can model coda wave envelopes depending on frequency, lapse time, and scale of the heterogeneities encountered during propagation (Sato et al., 2012). Coda waves can enter the diffusive regime at late lapse times, where coda attenuation is theoretically equal to absorption (Q c ≈ Q i - Shapiro et al., 2000).  demonstrated that this assumption is valid for an onset time of the coda of 80 s and epicentral distances between 0 and 90 km across the Pyrenees. Q c −1 strongly depends on the frequency in active tectonic regions (Sato et al., 2012). At low frequencies, high Q c −1 anomalies map surface geology and sedimentary basins in the Pyrenees , the Alps (Mayor et al., 2016), and Vrancea (Borleanu et al., 2017). At higher frequencies, the change in coda composition from surface to body waves can increase depth sensitivity (De Siena et al., 2016;Mayor et al., 2016;Gabrielli et al., 2020). This effect is dominant when measuring Q c −1 across volcanic structures, as documented along the Mount St. Helens volcano (US) where low-frequency waves allow to map the shallowest, most-heterogeneous volcanic structures, while deep feeding systems appear as high-Q c −1 anomalies at higher frequencies (De Siena et al., 2016;Gabrielli et al., 2020). Di Martino et al. (2022) recently obtained absorption maps using an active seismic experiment (~100 m 2 array with meter-scale station distance) at the Solfatara volcanic crater, in Southern Italy. This study demonstrates that absorption anomalies mark zones of high soil temperatures and CO 2 -bearing fluids migrating from an extended fault, whose features were reconstructed using peak delay maps.
High peak delays and high Q c −1 detect intrusions in mountain chains  and magmatic systems (De Siena et al., 2016). These parameters are reliable markers of fluid-induced fracturing leading to earthquakes at the regional scale (Borleanu et al., 2017). In fault networks like those ruptured by the AVN, Napolitano et al. (2020) found that high-scattering and highabsorption patterns obtained by increasing frequencies map the migration of the historical seismic events from the 16th century until the Pollino seismic sequence. This work supported the view that lateral fluid diffusion drove the sequence (De Matteis et al., 2021;Napolitano et al., 2021). Focal mechanism tomography highlighted the excess of pore pressure in the hypocentral area, which affected the seismicity distribution during the seismic sequence (De Matteis et al., 2021). In the Alps, at low frequency, high-absorption regions correspond to Cenozoic and Mesozoic sedimentary basins (Mayor et al., 2016). The Q c −1 map of the western Tibet shows variations along the Karakoram fault, attributed to the heterogeneous crust, velocity variations, and mantle-derived fluids (Biswas and Singh, 2020).
However, the temporal potential of joint scattering and absorption mapping in fault networks is still unexplored.
This study measures and maps peak delays and Q c −1 in the 2D space. Due to the consistent spatial extent of the earthquakes, we repeat the mapping in time using data recorded: The results offer insight into fluid migration and expansion processes leading and controlling the seismic sequence using peak delay and coda attenuation as frequency-dependent proxies of scattering loss and absorption.

GEOLOGICAL AND SEISMOLOGICAL BACKGROUND
The Central Apennine is governed by the opening of the backarc Tyrrhenian basin to the west and the eastern migration of the compressive front associated with the lithospheric plate's flexural retreat below the Italian peninsula (Malinverno & Ryan, 1986;Doglioni, 1991;Carminati et al., 2014). The western Tyrrhenian basin is marked by thin crust (<25 km), high heat (>200 mWm −2 ), and positive gravity anomalies. Most of the Quaternary volcanoes and geothermal areas (such as Larderello, Amiata, Campi Flegrei, and Vesuvius) are activated by the uprising of the mantle wedge between the Tyrrhenian Sea and the southwestern subducting Adriatic domain below the Apennines (Ventura et al., 2007). On the contrary, the Adriatic domain is characterized by low heat flow, negative gravity anomalies, and a crustal thickness of about 35 km (Carminati and Doglioni, 2012).
The Central Apennines were controlled by compressive tectonics (Galadini and Galli, 2000) up to Pleistocene times and later affected by an NE-SW extension. The Quaternary extension in the central sector of the chain is coeval with the eastward migration of the thrust fronts in the Adriatic foreland (Pizzi and Galadini, 2009). The NW-SE striking Central Apennine region is mainly composed of thrust sheets imbricated toward the Adriatic Sea (E-ENE direction) and creating a contact between the Meso-Cenozoic carbonate succession and the Miocenic flysch (Billi and Tiberti, 2009;Cosentino et al., 2010). The sequence of carbonate rocks and terrigenous units can be traced from the surface to an 8-10 km depth, moving from east to west (Billi and Tiberti (2009) and references therein). The Monti Sibillini thrust in the southern sector of the Umbrian Arc has an NNE-SSW trend. It is responsible for the overthrust of the Sabina Jurassic-Miocene pelagic sediments on the Lazio-Abruzzi carbonate platform (Billi and Tiberti, 2009). The Laga Formation (consisting of Messinian siliciclastic foredeep deposits) overlies both the Triassic-Miocene shallow water carbonates of the Lazio-Abruzzi platform and the Umbria-Marche pelagic successions ( Figure 1) (Cosentino et al., 2010;Scisciani et al., 2014;Di Bucci et al., 2021, and references therein).
Further east, the contact between the Laga Formation and the carbonate platform coincides with the E-W thrust of the Gran Sasso massif, bounded by the Monti Sibillini thrust to the west (Billi and Tiberti, 2009). Large intermountain basins as those of Amatrice, Norcia, and Castelluccio are filled by Plio-Quaternary continental sediments (Figure 1). They formed in response to the activity of Plio-Quaternary NW-SE striking normal faults dissecting the Apennine chain and partly reactivating the preexisting, older thrusts Scognamiglio et al., 2018;Improta et al., 2019, and references therein).
Seismicity in the Apennines concentrates along the chain axis and manifests itself through swarms and sequences with events of magnitudes up to 7 and depths shallower than 10-15 km (Chiarabba et al., 2005;Chiaraluce et al., 2004Chiaraluce et al., , 2017. The mainshocks have dip-slip, normal focal mechanisms consistent with ruptures along NW-SE striking fault planes moving in response to an NE-SW extension (Montone and Mariucci, 2016;Buttinelli et al., 2018;Scognamiglio et al., 2018). Although some mainshocks nucleate in the crystalline basement underlying the chain at about 10 km depth, there are records of earthquakes in the overlying carbonates and sedimentary sequences (Chiaraluce et al., 2017;Buttinelli et al., 2018). Although the seismicity of Apennines is related to the extensional tectonics affecting the chain from Plio-Pleistocene times, the release of deep-derived, CO 2 -rich fluids has been proposed as a possible concurrent triggering mechanism of earthquakes (Miller et al., 2004;Di Luccio et al., 2010;Malagnini et al., 2012;Chiodini et al., 2004Chiodini et al., , 2020. Chiodini et al. (2020) showed a record of 10 years (2009-2018) of tectonic CO 2 release temporally correlated with the number and energy of earthquakes in the Central Apennines (Italy), where devastating historical earthquakes (such as the 1461 event of L'Aquila, the 1703 event of Norcia-Montereale-L'Aquila, and the 1915 M w 7.0 event of Avezzano) occurred. Chiaraluce et al. (2007) and Collettini et al. (2008) reported over-pressurized fluids in the deep wells of San Donato and Santo Stefano in the Apennines. Miller et al. (2004) proposed that the 1997 Colfiorito aftershocks were caused by a co-seismic release of trapped high-pressure CO 2 . Chiarabba et al. (2009b) observed high V p /V s during the foreshocks and aftershocks of the L'Aquila 2009 sequence (M w 6.1), suggesting high pore pressure and fluid storage zones at depth. Malagnini et al. (2012) andDi Luccio et al. (2010) showed the diffusion of fluids and pore fluid pressure variations during the L'Aquila earthquake through the spatial and temporal evolution of seismicity.
In the Apennines, earthquakes with M > 6.0 often rupture the surface, so fluids circulating in the shallower crust can rise and interact with the fault zones (Amoruso et al., 2011;Doglioni et al., 2014). The Central Apennines host large aquifers in the carbonate formations, such as the Gran Sasso (Amoruso et al., 2011) and Nuria-Velino-Giano (Devoti et al., 2018) aquifers. The Gran Sasso aquifer registered short-and mid-term changes in groundwater hydrodynamics (Amoruso et al., 2011) after the 2009 L'Aquila seismic event. The Gran Sasso aquifer extends 700 km 2 , bounded by the Laga Formation to the north and east and by low-permeability alluvial deposits to the west and south. The perennial groundwater reserves in this aquifer are estimated in the order of 10 10 m 3 , with a mean thickness of 1 km (Amoruso et al. (2011) and references therein). The 2016 Amatrice-Norcia seismic sequence generated an uplift of the water level near Gran Sasso (~1.8 m; Devoti et al., 2018) and 100 km away from the mainshock area (up to 80 cm; Barberio et al., 2017). The Mounts Nuria-Velino-Giano hydrological complex is located~40 km southwest of the Gran Sasso aquifer and has an outcropping area of about 1,000 km 2 (Boni et al., 1987;Chiodini et al., 2011). It is bordered by low-permeability deposits (Cretaceous-Pliocene flysch formations) and the Mt. Sibillini thrust (Boni et al., 1987;Chiodini et al., 2011). Chiodini et al. (2011) performed a hydrogeochemical study of the Gran Sasso and Nuria-Velino-Giano aquifers to investigate the existence of a deep CO 2 source in the epicentral area of the 2009 L'Aquila seismic sequence. Both aquifers showed an increased influx of fluids rich in deep CO 2 before and during the seismic sequence. Thus, deep high-pressure gas traps played a role in triggering the L'Aquila seismic sequence, as previously modeled for the Colfiorito seismic event (Miller et al., 2004). Tomography images suggest that such traps occur at the base of the Apennine seismogenic layer, i.e., at the boundary between the upper and the lower crust at 10-15 km depth (Chiarabba et al., 2020a).
The AVN sector of the Central Apennines is~30 km north of the L'Aquila area and~50 km south of the hypocenter of 1997 M w 6.0 Umbria-Marche seismic sequence. The seismicity of this sector is associated with the Mt. Vettore normal fault systems and the Mt. Gorzano fault (Monti della Laga fault system) Brozzetti et al., 2019;Carminati et al., 2020). The AVN ruptures fill the seismic gap between the 1997-98 Colfiorito  and 2009 L'Aquila seismic sequences (Chiaraluce et al., 2011).

Seismic Dataset
In the present study, we used a large set of seismograms merging weak-and strong-motion data from the AVN using earthquakes with magnitudes between 2.8 and 6.5. We selected earthquakes having a maximum depth of 20 km, keeping source-station distance within 100 km.
Strong ground motion data were registered by the accelerometric stations of the Italian Strong Motion Network (RAN). Broadband weak-motion seismological stations are part of the Digital Seismic Network run by Istituto Nazionale di Geofisica e Vulcanologia. The weakand strong-motion accelerograms were downloaded from the ITalian ACcelerometric Archive (ITACA) website, the European Strong Motion (ESM) database, and the European Integrated Data Archive (EIDA) repository. The strong ground motion network is equipped with a threecomponent Kinemetrics EpiSensor (FBA-3200 Hz) with a full-scale range of 1 or 2 g, combined with ETNA 18 bits or K2-Makalu 24 bits digitizers. The weak-motion seismograms are corrected for instrument response. All the stations that registered only one seismic event were removed. The waveforms with P-wave travel times higher than 35 s were discarded from the database to constrain their propagation within the crust. Waveforms with spikes, telemetry gaps, and wave arrivals in the coda of the horizontal components of ground motion were manually removed. The dataset comprises seismograms with a signal-to-noise ratio (SNR) of the selected coda window higher than three. To estimate the noise level, a 3 s time window is considered before the P-wave arrival. The P-and S-wave arrivals were picked manually on each seismogram of the final dataset.
The final dataset comprises~3,200 waveforms for the presequence ( Figure 2A) recorded at 40 seismic stations and 13,600 recorded at 154 seismic stations ( Figure 2B) for the AVN. The number of waveforms for the Q c −1 measurements decreased to~1,200 recorded at 39 stations for the pre-sequence and~5,800 recorded at~141 stations for the AVN due to the signal-to-noise selection. The number of stations working during the sequence is almost four times larger than the one recorded during the pre-sequence due to a temporary seismic network deployment just after the Amatrice mainshock of August 24. We selected similar locations of events and source-station pairs considering the expected resolution of our final maps to carry out an appropriate comparison of the analyses of the two datasets in different periods. In the following sections, we investigated and demonstrated that the location of source-station pairs or difference in the network configuration densified during the sequence does not significantly impact the spatial pattern of scattering and absorption in the studied area.
Our study covers an area of 200 km × 220 km (lon: 11.85°-14.25°; lat: 41.7°-43.7°), and it is divided into regularly spaced nodes. The seismograms were filtered with a band-pass Butterworth filter (fourth order), forward and backward, in four frequency bands (1-2 Hz, 2-4 Hz, 4-8 Hz, and 8-16 Hz) centered at a frequency (f c ) equal to 1.5, 3, 6, and 12 Hz. A Hilbert transform has been applied to compute the envelopes, which were then smoothed with a moving window of duration eight times the central frequency's inverse.
The open-access code MuRAT, previously applied in volcanic (De Siena et al., 2016Gabrielli et al., 2020;Di Martino et al., 2022) and tectonic settings (Borleanu et al., 2017;Napolitano et al., 2020), is used to perform all the analyses. The code is applied to datasets spanning five time frames: the presequence (2013-2016), the AVN, and the three sub-sequences after the mainshocks (from August 24, 2016, to January 18, 2017). Peak delay and Q c −1 parameters mapped in space using the MuRAT computer codes are then mapped using the Generic Mapping Tool, GMT (Wessel and Smith, 1998), and interpreted separately in their parameter space in the following sections.

Peak Delay Measurements and Mapping
Seismic waves broaden in heterogeneous media due to multiple forward scattering, with scattering attenuation affecting the earlier parts of seismic envelopes (Zhang et al., 2021). The delay between the S-wave onset and the maximum of the envelope, which increases with the source-receiver distance, is thus a proven marker of scattering attenuation (Takahashi et al., 2007;. The dependence of the peak delay times (t PD (r), in seconds) on hypocentral distances (R Hypo , in km) is expressed by where A(f) and B(f) are the coefficients of the resulting linear fit (Gusev and Abubakirov, 1999). We selected seismograms with an R Hypo between 20 and 100 km and calculated empirical peak delay times against hypocentral distances in logarithmic scale for four frequency bands since the envelope broadening for waveforms close to the source (R Hypo < 20 km) may be caused by the source duration and prevails over scattering effects (Takahashi et al., 2007). Measured peak delay times over the E-W horizontal component of the seismograms (black circles) and their dependence on hypocentral distances are presented for the frequency band 1-2 Hz in Figure 3 for all sequences analyzed in our study. While peak delays show significant scattering, they exhibit a linear trend and increase with distance. Red lines in Figure 3 mark the linear regression lines of peak delay times against hypocentral distance in the frequency band 1-2 Hz, which can be characterized as the theoretical peak delay manner. The difference between the measured peak delay time of the ith waveform (t PD i (f)) and the theoretical peak delay at the corresponding hypocentral distance (Eq. (1)) gives the amount of scattering accumulated along the ray path: Positive values of Δ log 10 t(f) represent high-scattering zones and heterogeneous portions of the crust. In contrast, negative values mark low-scattering zones interpreted as rigid and compact crustal areas. Figure 4 shows some seismic envelopes observed at two stations located at the same distance from events, 35 km. The envelopes from the two selected events are shown for the 1-2 Hz frequency band. It is acknowledged that the envelope FIGURE 3 | Log-log plot of peak delay (s) as a function of the hypocentral distances (km) for the pre-sequence, 2016-2017 sequence, and single sequences at 1.5 Hz. The red solid line represents the best linear fit for the peak delay measurement (black circles). A(f) and B(f) are the coefficients of the best linear fit.
Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 909698 on ray path R1 reveals larger peak delay times than those of ray path R2. This statement strongly implies the presence of strong and weak heterogeneities along with R1 and R2, respectively. Peak delays thus indicate the presence of scattering in the rock volumes encompassed by the propagation of the primary shear waves. The size of such volumes depends on the wavelength and on the amount of scattering and should be estimated by appropriate sensitivity kernels. The assumption of a source-receiver sensitivity "on ray" produces a large underestimation of the size of these volumes, which in reality are larger than those associated with the primary waves, sometimes named "banana-doughnut" (Nolet, 2008). This on the one hand makes the space sensitivity of peak delay imaging lower than that of travel time imaging but on the other hand should reduce the bias introduced by changing the input data in the ordinary travel time and attenuation tomography.
We are thus confident that peak delay mapping produces estimates of the true scattering attenuation patterns which should be more stable and robust than those obtained using ordinary methods.
To quantitatively portray the observed spatial variation in peak delay times at stations, we mapped the logarithmic deviation of S-wave peak delay time, assuming source-receiver sensitivity on rays and using a standard regionalization approach. We divided the study region into many small blocks of different sizes horizontally and in depth. So, we considered the spatial variation of medium heterogeneities by assigning the averaged value of Δ log 10 t(f) to each block. The ray coverage is dense in the central sector of the map; however, peak delay values spike in regions of low ray coverage. To avoid these trade-offs in areas surrounding the seismogenic zone, we removed all nodes crossed by less than ten rays . We performed several tests on the peak delay parameter and its sensitivity and variability to the location of source-station pairs or the network configurations, number of ray paths, and block size to select appropriate parameters and support our results' resolutions in space and time. First, we divided the study area into 15 × 15 and 30 × 30 blocks, of which the size is 0.12°× 0.12°and 0.06°× 0.06°horizontally and 20 km in depth. We then changed the minimum number of ray paths that travel each block from at least 10 to 20 and more ray paths. We observed that both the grid size and the number of ray paths selected seem to have a negligible impact on its spatial pattern in the study region. The detailed presentations of these results can be found in the supplementary material in Supplementary Figure S1.
In addition, to investigate the peak delay parameter variation in depth, we split the Norcia dataset into two subsets: one between 0 and 5 km depth (~4,000 waveforms) and a second between 6 and 18 km depth (~3,500 waveforms), and compare the peak delay spatial variations from the two different datasets to investigate potential differences in the depth extent of the seismicity in different phases (Supplementary Figure S2). We stress that these events have been relocalized by Waldhauser et al. (2021). The results show minimal differences in peak delay mapping between the two datasets.
Finally, the station variability between the events registered during the pre-sequence and sequence is also examined in this section. To do so, we selected the same dataset of Norcia and reorganized the events recorded by the 40 stations of the presequence and the 154 stations of the entire sequence. We compared the spatial distribution of peak delay determined from different databases, decreasing 20% of the primarily selected events recorded during the sequence. However, Supplementary Figure S3 reveals how the peculiarities of peak delay are nearly identical in the two maps and qualify the robustness of our results and comparisons accomplished in space and time.
We also analyzed peak delay distributions during the Norcia sequence for the three components of motions at all the frequencies to study possible alterations that the anisotropy may cause. In Supplementary Figure S4, we compare the distribution of logarithmic deviation of S-wave peak delay time Δlog 10 t(f) measured on different components of seismograms. While the results obtained using the E-W and vertical components are almost identical, the N-S component shows opposite anomalies, corresponding to high scattering in the northernmost area of the seismogenic zone of Amatrice-Visso-Norcia. The slight contrast in the anomaly pattern corresponds to the Monti Sibillini thrust and the Syncline Camerate valley. This imprint observed at low frequencies can be explained by guided waves trapped between the Monti Sibillini thrust and the Syncline Camerate valley.
We prepared our final maps for detailed and better visualization of peak delay variation in the area by dividing the study area into 30 × 30 regularly spaced nodes and considering at least ten ray paths that travel each block. Supplementary Figure S5 shows the hit count maps where the blocks crossed by ten or more ray paths are colored in black for the pre-sequence, 2016-2017 sequence, and three sub-sequence datasets. Since our study focuses primarily on 2D spatial variation of S-waves' peak delay time, we considered all events to have depths between 0 and 20 km, entirely within the crust above the Moho, with no depth dependence for the resultant peak delay mapping. The results obtained based on the preferred parameters are given and discussed in Section 4.1 and Section 5.1 in detail.

Coda Attenuation Measurements and Mapping
The latter part of a seismogram is primarily sensitive to seismic absorption (Zhang et al., 2021). Aki and Chouet (1975) showed that the power spectral energy density (E(t,f)) is a function of the lapse time from the earthquake's origin time (t): where S(f) includes both source and site terms and Q c −1 is the frequency-dependent inverse coda quality factor. α is equal to 3/2 if propagation is constrained within a single layer characterized by an anisotropic multiple scattering regime (Paasschens, 1997;. Q c −1 is computed by linearizing Eq. 3. The linearization is computationally faster than non-linear gridsearch algorithms (e.g., Napolitano et al., 2020). The two approaches provide equivalent results when the signal-to-noise ratio is higher than three and for coda windows longer than 10 s (Sketsiou et al. (2020) and references therein).
Taking the logarithm of Eq. 3 and rearranging the terms, we can measure Q c −1 using a straight-line fitting: The Q c −1 analysis window starts at twice the arrival time of the S-wave arrival (tS) and lasts 20 s (length of the coda window), meaning that the central coda lapse time is on average higher than 3tS. To assess if this lapse time is sufficient to consider the waveforms as late coda, we determined the mean free time (MFT) considering the values of B 0 and Le −1 previously calculated by Akinci et al. (2020) through the multi-lapse time window analysis (MLTWA). Having the scattering coefficient η s = B 0 /Le, the mean free path (MFP) = 1/η s , and MFT = MTP/V s , we obtained a value of MFT = 9.8 s for 1.5 Hz and MFT = 26 s for 6.0 Hz. In the calculation, we assumed an average S-wave velocity in the crust of V s = 3.5 km s −1 . The average S-wave arrival time of the dataset is 39 s, with a minimum of 8 s and a maximum of 65 s: a coda window starting at 2tS thus contains primarily diffuse waves, allowing us to interpret Q c as a marker of absorption. The MFT can be calculated with more accuracy using the code Qopen (Eulenfeld and Wegler, 2016), which separates intrinsic and scattering attenuation for each event and station by envelope inversion. However, the comparison between Qopen and the MLTWA made by van Laaten et al. (2021) shows that both methods give the same values within their error bars. Figure 5 shows source-station measurements (black dots) and fitted (red line) Q c −1 obtained using the AVN dataset at two In addition, in this section, we presented results from several potential sources of variability included in the data, such as the difference in the network configuration, which was densified during the sequence. We used the same dataset of Norcia applying the identical criteria adopted previously for the peak delay parameter: considering events recorded by 1) the 40 stations of the pre-sequence and 2) the 154 stations of the entire sequence in our analysis. This process caused a decrease of waveforms from 3,200 to 2,300. However, Supplementary Figure S10 shows how the anomalies are essentially the same, allowing us to continue comparing the two time periods without excluding the temporary stations of the AVN dataset.
The depth variability of Q c −1 is also given in Supplementary Figure S11 from earthquakes distributed primarily across the seismogenic area, between 42.6°and 43°lat, for both the depth subdivisions. The seismogenic areas are dominated by a high attenuation anomaly, which appears more concentrated at higher depths rather than in the surficial dataset (Supplementary Figure  S11). This behavior can be ascribed to a more fractured shallow crust and a more compact and less fractured deep crust. Nevertheless, high coda attenuation is still focused in the central area independently of depth, so we excluded a dependence on the earthquake depth variations in this study.

RESULTS
The results are shown between 12.5°E-13.7°E and 42.15°N-43.48°N for the peak delay spatial variation and between 12.7°E-13.52°E and 42.37°N-43.48°N for the Q c −1 spatial distribution, as a zoom of the earthquake-damaged and the fractured area in Figures 2A,B. Areas of high resolution for peak delays and Q c −1 are defined using hit counts (Supplementary Figure S5) and the results of the checkerboard tests (Supplementary Figures S8, S9), respectively. We presented the individual spatial variations of peak delay and Q c −1 using different sets of maps. The final results and interpretations focus on how the scattering and absorption parameters vary in time and space in the study area. Figure 6 shows the spatial variations of the peak delay parameter (scattering losses, as defined in Eq. 2) obtained at 1.5 Hz for the pre-sequence, the AVN, and the three sub-sequences. The final peak delay maps are frequency-dependent and assume no dependence on depth, thus mapping heterogeneities of Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 909698 different dimensions. The earthquake dataset is entirely within the crust, between 0 and 20 km depth, so above the Moho (~35-40 km depth, Piana Agostinetti, and Amato, 2009), and the wave propagation is constrained above the primary horizontal interface. The maps show the absolute value of the frequencydependent peak delay, Δ log 10 t(f): low scattering and high scattering are the values negative and positive around zero, which is the average logarithmic peak delay of the area. The results at higher frequencies (3, 6, and 12 Hz) are shown in Supplementary Figures S12, S13, S14. Areas of small Δ log 10 t(f) values conform to weak inhomogeneities; lowscattering attenuation is marked by cold (blue) colors, while large Δ log 10 t(f) values, high-scattering attenuation, are hot (red) colors. The blocks crossed by more than or equal to 10 ray paths are shown in Supplementary Figure S5.  Figure 6). To the east of the carbonate mountain ridge, a positive anomaly is associated with Pliocene synorogenic formations on top of the carbonate thrust units Bigi et al., 2011). We observed variations in scattering attenuation over different periods. During the pre-sequence, the location of the Colfiorito earthquake of 1997 in the northernmost portion of the maps is predominantly low scattering. Highscattering attenuation characterizes the area during the Amatrice and Visso sequences, and the area becomes progressively low scattering during the Norcia sequence. This area is delimited to the south parallel to the Monti Sibillini thrust (MST) (Figures 1, 6). Before the Amatrice earthquake, the area intersected by the AVN was marked by low-scattering values. During the Amatrice sequence, the southeastern part of the epicentral area becomes a low-scattering region, while scattering increases northwest. In general, scattering losses increase and expand toward the epicentral area during the AVN. This behavior is evident when comparing the presequence with the entire AVN ( Figure 6A). The lowscattering anomalies progressively reach the southwestern patch of the Amatrice fault rupture area, spreading across the dense fault network west of Norcia during the Visso and Norcia sequences ( Figure 6B). At higher frequencies, a high-scattering anomaly appears west of the mountain range (Supplementary Figures S12, S13, S14).

Average and Spatially Varying Q c −1
We report the frequency-dependent average Q c −1 and their relative standard deviation for the pre-sequence and the AVN datasets in Table 1. These values are obtained from the best fit line of Eq. 4. Q c −1 decreases with frequency, as expected in tectonically active regions (de Lorenzo et al., 2013;Akinci et al., 2020). Assuming a frequencydependence given by Q c −1 = Q 0     Figure S15). High-absorption anomalies are widespread before the Amatrice mainshock and focus on the fault zones during the AVN. The L'Aquila-Gran Sasso sector, GST, shows a drastic decrease in attenuation during the AVN. High absorption concentrates on the lower fault plane during the pre-sequence, shifting further north after the Amatrice mainshock (red rectangles in Figure 7). High absorption marks the dense fault network south of Norcia and the fault planes after the Visso mainshock. After the Norcia mainshock, the NNW-SSE-trending high-absorption anomaly expands across the Central Apennine range, bounded north and south by the Monti Sibillini thrust (MST) and Gran Sasso thrust (GST), respectively. The dense fault network south of Norcia shows a substantial absorption decrease at the end of the AVN compared to the Amatrice and Visso periods. We identify a similar trend at high frequencies (3, 6, and 12 Hz), with a more focused high anomaly during each seismic sequence (Supplementary Figures S15, S16, S17).

Parameter
To demonstrate the robustness of our interpretation of temporal evaluation of absorption at a smaller scale and to examine the distribution of the anomalies as a function of time, we divided the Norcia dataset into two subsets: the first subset between October 30 and November 7, 2016, and the second between November 7, 2016, and January 18, 2017. Each dataset contains 245 events. As demonstrated in Figure 8, the distribution of the earthquakes between the two datasets is similar, apart from the southernmost area, more covered by the events during the second subsequence (below 42.6°lat). The locations of the absorption pattern anomalies are different between the first and the second subset of the dataset and change with time, even during the Norcia sequence. This analysis reveals that coda attenuation evolves in the seismogenic zones independently of earthquake distribution, allowing us to interpret the maps in terms of temporal processes. This is expected as coda waves become progressively independent of source-station locations at late lapse times.
Other than mapping inverted Q c −1 value deviation in time and space, we also compared the averaged Q c −1 at two seismic stations (GUMA and LNSS) located in the northern and southern sections of the ruptured fault area, respectively, where substantial changes in absorption were found between January 2013 and 2017. In Figure 9, we plot averaged Q c −1 values calculated at two stations together with the magnitude of the seismic events of our dataset during the same period and the deeply derived carbon (C deep ) at the Canetra spring and the dissolved CO 2 in the groundwaters of the Velino aquifer (FCO 2 ), both measured by Chiodini et al. (2020). Figure 9A denotes an apparent change of absorption before and after the mainshock of Amatrice, while Figure 9B better highlights these changes during the sequence period. During the presequence (from January 2013 to August 2016), absorption fluctuates between 0.0050 and 0.0070 at both stations ( Figure 9A). However, it varies more rapidly and significantly (between 0.0035 and 0.0080) during the sequence period (approximately 5 months) ( Figure 9B). Although the most abrupt changes in absorption correspond to the three mainshocks (quick drops and peaks), the general trend is different at two stations, decreasing or increasing over the whole time of the sequence. Notably, the absorption always decreases quickly a few days after the Amatrice mainshock at station LNSS, located in the southwestern part of Amatrice; on the contrary, it increases at station GUMA, located in the northeastern part of the faulting area, closer to Visso. There is also a quick drop in Q c −1 after the Amatrice earthquake, followed by a rapid increase till October 26 and 30 at station GUMA. After the Norcia earthquake, a reduction in absorption is evident at both stations.

DISCUSSIONS
In the following, we discussed peak delay and Q c −1 variations in space and time separately (Figure 6 and Figure 7). These variations are investigated as a function of time and space for the pre-seismic and seismic periods. Spatial differences between the pre-seismic and the seismic phases occur because of different processes, which include fluid flow, crack opening/closure and crack density, pore pressure variations (compaction/dilatation), lithological contrasts, and existing fault networks. These processes produce different attenuation mechanisms.

Peak Delay (Scattering) Patterns in Space and Time
At all frequencies and for all sequences, a low-scattering NNE-SSW directed anomaly coincides with the Monti Sibillini thrust, including the flysch of Laga Formation. Thus, the thrust appears as a low-scattering barrier for northern propagation of fluidinduced seismicity (Improta et al., 2019). High-scattering losses correspond to the carbonatic Gran Sasso massif and the L'Aquila basin (Chiaraluce et al., 2011), where the 2009 L'Aquila seismic sequence nucleated. These losses mark high velocity and density contrasts, high density of cracks and fractures, and basins in the upper part of the crust (Borleanu et al., 2017). The correspondence between high scattering and carbonate rocks fractured by historical sequences is analogous to that observed across the Lauria Mountains in Southern Italy . Considering the correlation with the Gran Sasso and the Monti Sibillini thrust (GST and MST; Figure 6), peak delays appear primarily sensitive to the existing structural elements as the pre-existing thrust and active faults De Siena et al., 2016). This inference is supported by the spatial correlation with the velocity model of the Central Apennines of Chiarabba et al. (2020a). The authors show lateral heterogeneity across the Monti Sibillini thrust, with high V P and V P /V S anomalies consistent with the highscattering Lazio-Abruzzi domain and low V P and V P /V S in the northern portion corresponding to the lowscattering Umbria-Marche domain. The V p peculiarities characterize a primary discontinuity between the Monti Sibillini and Laga thrust units, with the uplifted sedimentary deposits forming along distinguishable thrust faults (Centamore et al., 2009).
Peak delay patterns change less than Q c −1 patterns in time and frequency ( Figure 6 and Supplementary Figures S12, S13, S14). However, while high-scattering anomalies appear diffuse before the AVN at 1.5 Hz, they focus on a narrow NW-SE striking band just south of the fault planes during the sequence. Moving from the Amatrice to the Norcia sequences, the high-scattering anomaly progressively expands from the southern fault plane to comprise the dense fault network located immediately to the W-SW of the modeled faults (red boxes in Figure 6). We interpret this change in scattering attenuation as evidence of micro-fracturing processes within the fault network, lacking evidence of significant lithological changes. This inference is supported by the higher-frequency peak delay maps (Supplementary Figures S13, S14) and by the high-anisotropy values measured by Pastori et al. (2019) in the same area and interpreted as heavily fractured, stressed, and overpressurized by channeled fluids during the AVN.

Coda Q (Absorption) Patterns in Space and Time
The sparse distribution of the high-absorption anomalies during the pre-sequence becomes a continuous NW-SE striking pattern focused on the seismogenic zone during the AVN ( Figure 7A). Spatial variations of high-absorption patterns are the manifestation of increased fluid content, with fluids expanding across the seismogenic zone with each sequence ( Figure 7B). Fluids migrate across the faults at seismogenic depths (8-12 km), where the main events of the AVN have occurred. Absorption is higher in the southern part of the ruptured area during the Amatrice earthquake sequence; then, the high-absorption pattern moves north during the Visso sequence and extends along the fractured zone during the Norcia sequence. The highabsorption pattern thus points to northern migrations of fluids and across tectonic structures associated with fracturing up to the Norcia epicenter. At the end of the sequence, absorption has decreased on average, especially in the southern section of the fault zone where the Amatrice earthquake occurred. The images suggest that the fluids had abandoned the area by the end of the seismic sequences.
In fluid-filled fault networks, high absorption marks regions with the highest concentrations of CO 2 -bearing fluids (Di Martino et al., 2022). Our absorption patterns evolve during the sequence as expected for high-absorption deep CO 2 -bearing fluids previously trapped in the fault systems (Miller et al., 2004). The development of a high-pressure front allows fluids to propagate through the seismogenic zone and create aftershocks. The absorption patterns show the evolution of this front from the location of the Amatrice earthquake toward the north and south when fluids diffuse across the fault planes (Figure 7). Northern migration happened primarily during the Visso sequence when Monti Sibillini acted as a barrier (Improta et al., 2019). The Norcia sequence is the manifestation of the expansion of the over-pressurized trapped fluids across the fractured region (Pastori et al., 2019), marked by an expanded absorption anomaly.
The structural models confirm that fluid pathways follow the maximum horizontal stress direction (Sibson, 2000), which in the Central Apennines is NW-SE-oriented and parallel to the strike of the active normal faults (Di Luccio et al., 2010). The rock formations on the western side of the fault system are heavily fractured, and apparently, they channeled and trapped fluids. High-absorption anomalies focus primarily west of the activated fault system during the AVN (Figure 7). They mark the extension of the highly fractured and fluid-filled zones and the area of maximum displacement recognized by Brozzetti et al. (2019) for the Monte Vettore fault. This pattern is even more evident in Figures 8A,B, where after less than 1 month from the Norcia mainshock, the high attenuation anomaly is all focused east of the sequence fault plane ( Figure 8B).
Our observations are consistent with the results from geochemical and geophysical studies that infer the critical role of CO 2 -rich deep fluids in the seismogenesis and evolution of seismic sequences across the Apennines (Chiodini et al., 2004(Chiodini et al., , 2020Miller et al., 2004;Di Luccio et al., 2010;Terakawa et al., 2010;Malagnini et al., 2012;Barberio et al., 2017;Chiarabba et al., 2020b). These studies highlight a relation between CO 2 release and seismicity rates and locations (Miller et al., 2004;Chiodini et al., 2020). Such CO 2 is supposedly released mainly from the mantle (Chiodini et al., 2004(Chiodini et al., , 2020Frezzotti et al., 2009) and accumulates at the lower and upper crust boundaries. Here, overpressurized reservoirs develop, sealed by 10-12 km deep evaporitic levels (Chiodini et al., 2004(Chiodini et al., , 2020, and by the multilayered carbonate succession in the inherited thrust (Chiarabba et al., 2020b;Carminati et al., 2020). Chiarabba et al. (2020b) provided temporal evidence of a reduction of the V p in the vicinity of the mainshock nucleation. They associated this variation with increasing fluid pressure, lubrication, and weakening of the fault, as confirmed by laboratory experiments (Di Toro et al., 2011). Soldati et al. (2019) investigated the temporal and spatial variation of the crustal velocity during the seismic sequence from noise cross-correlation. They observed a sudden drop in the seismic velocity at each mainshock, followed by a quick recovery. After the sequence, the co-seismic velocity gradually returned to the presequence level, with a 75% perturbation recovered by the end of 2017. In space and time, this velocity drop was interpreted as the effect of the fast circulation of fluids, more abundant in the L'Aquila area, where the carbonate platform is located. The velocity drop distribution fluctuates with time and space, with the highest variation along the NE direction, corresponding to the sandstones of the foredeep deposits, highly hydrated and porous. These deposits increase pore pressure during the seismic sequence, decreasing the effective relative velocity.
Indeed, an increase in the concentration of mineralized endogenic fluids corresponds to seismic sequences (Barberio et al., 2017). The source of these fluids has been identified in rocks located at a depth of about 10 km (Figure 1). The influx of deep CO 2 could favor their formation and ascent to the surface. According to Di Luccio et al. (2010), these deep fluids mainly migrate along the NW-SE striking faults and fractures of the Central Apennines. However, our results do not exclude those fluids moving in the shallower portion of the crust, especially where the permeability is higher, as within carbonates and basin filling deposits. Regardless, the maximum recorded water level variation in the Apennine aquifers during the AVN (e.g., Gran Sasso aquifer) is only 1.8 m (Devoti et al., 2018).
The timelines in Figure 9 show the time evolution of averaged Q c −1 variation next to the main fault(s) or at those locations (stations GUMA to the north and LNSS to the south) where the most striking changes are observed. Figure 9 helps evaluate the critical variations over the pre-sequence time and remark the details during the sequence period. In general, we observed a smooth and steady absorption trend over the three and a half years of pre-sequence. This variation is followed by a rapid rise and drop of the absorption starting from the Amatrice mainshock, which increased and diminished after a few days. Significant absorption changes are mainly detected before and after the three mainshocks of the sequence.
A similar pattern was observed by Chiodini et al. (2020) in the temporal evolution of deep CO 2 degassing. Furthermore, Soldati et al. (2019) investigated the velocity perturbations from crosscorrelation of ambient seismic noise associated with the Central Italy seismic sequence of 2016-2017. According to the individual of the three main shocks of the sequence, they observed an abrupt reduction of seismic velocity, followed by a quick recovery in the short term. They associated these fluctuations with fluid migrations and pore pressure changes and their effect on crack dynamics (opening-closure). The rapid change of absorption between Visso (October 26) and Norcia (October 30) is also identified in the 4D velocity tomography of Chiarabba et al. (2020b), where a decrease in V p /V s occurred near the Norcia hypocenter.
Our findings on the absorption behavior in time and space are consistent and comparable to those documented from independent studies using multi-parametric data in the same region. The decrease in absorption at the LNSS station and the increase at GUMA right after the Amatrice earthquake may emulate the fluid movements from south to north over the ruptured area. After the Amatrice earthquake, fluids moved to the north and triggered Visso, where we observed a significant increase in absorption at GUMA while a significant decrease in absorption at LNSS where the fluids have abandoned the area moving to the north.

CONCLUSION
The main results of this study are the following: 1) Peak delay analysis delineates and characterizes the structural (e.g., Monti Sibillini thrust) and lithological (e.g., fractured carbonate rocks) domains of the Central Apennines over all the studied frequency bands. High-scattering patterns highlight fracturing occurring during the AVN, with highscattering anomalies focusing south and west of the seismogenic zone. 2) Coda-Q mapping marks the concentration of fluids across the seismogenic zone during the sequence. It detects a progression in time associated with fluid expansion from the Amatrice epicenter, mapping the migrations of deep CO 2 -rich fluids across the fractures and the fault networks of the Central Apennines.
These results point to a common feature for earthquake sequences across the Central Apennines: the deep migration and expansion of CO 2 -bearing fluids along the NW-SE elongated seismogenic zone. The front and area of expansion of these processes can be detected by mapping different seismic attenuation mechanisms in time and space. Knowing their variations also profoundly impacts assessing seismic hazards from earthquake triggering to earthquake ground motion to mitigate risks.
Mapping scattering and absorption in the 3D space through seismic sequences could give us unprecedented sensitivity to structures and processes leading to mainshocks and