The Attenuation and Scattering Signature of Fluid Reservoirs and Tectonic Interactions in the Central‐Southern Apennines (Italy)

Despite the high detection level of the Italian seismic network and the risk associated with its fault networks, Central‐Southern Italy has no unique geophysical model of the crust able to illuminate its complex tectonics. Here, we obtain seismic attenuation and scattering tomography models of this area; both reveal high attenuation and scattering anomalies characterizing the entire Apenninic Chain and related to its East‐ and West‐dipping extensional Quaternary tectonic alignments. Fault‐associated fractured zones become preferential ways for circulating and degassing high‐attenuation CO2‐bearing fluids. A previously undetected fluid source area is a high‐attenuation volume below the Matese complex, while a similar smaller anomaly supports a fluid source near L'Aquila. The most prominent low attenuation and scattering volumes reveal a locked aseismic zone corresponding to the Fucino‐Morrone‐Porrara fault systems, representing a zone of significant seismic hazard.

In this context, shallow seismic imaging can reveal crustal geometries at the first tens of kilometers, connecting the lithospheric context with the surface or sub-surface geology.However, no seismic crustal image of the Central-Southern Apennines transition zone is currently available.Published tomographic models only reach the area marginally (Chiarabba et al., 2010;Improta et al., 2014;Zhao et al., 2016).This is primarily due to the low rates of seismic activity (Bagh et al., 2007;Frepoli et al., 2017;Romano et al., 2013;Trionfera et al., 2019).Nevertheless, the increased detection capability of the Italian permanent seismic network during the last two decades provides sufficient data for obtaining a regional-scale tomographic crustal model of the area.
Attenuation tomography techniques use seismic amplitudes from earthquake recordings to measure and model the energy lost by P-and S-waves while propagating through space (Haberland & Rietbrock, 2001;Lees & Lindley, 1994).Especially effective when targeting fluid and melt propagation at the lithospheric scale, the technique has been remarkably efficient in imaging melt sills and fluid reservoirs in volcanoes (De Siena, Thomas, & Aster, 2014;De Siena, Thomas, Waite, et al., 2014;De Siena et al., 2017;Prudencio, De Siena, et al., 2015;Prudencio, Ibáñez, et al., 2015).Its most recent applications to fault networks and mountain ranges show increased sensitivity of the technique to fluid-driven processes and strain accumulation (Amoroso et al., 2017;Sketsiou et al., 2021), as expected from laboratory experiment and modeling results (Di Martino et al., 2022;Tisato & Quintal, 2014).
Fracture networks attenuate body waves primarily by scattering (King et al., 2022;Sato et al., 2012), so much so that scattering tomography is increasingly used to map structural barriers, like thrusts and significant faults, across mountain ranges and volcanoes (Gabrielli et al., 2023;Napolitano et al., 2020;Reiss et al., 2022).Coda waves are the signature of such scattering processes (Sato et al., 2012).With increasing heterogeneity, codas build up, broadening waveforms, and delaying the peak of seismic envelopes (Saito, 2002).Peak delays increase abruptly at stations near shear zones due to trapped and resonant waves, providing an excellent marker of these features (King et al., 2022(King et al., , 2023)).Late lapse-time coda wave energies can also normalize direct-wave energies (coda normalization), removing source and site trade-offs from path-dependent attenuation measurements (Sato et al., 2012).The coda-normalization method has thus become a standard in attenuation tomography in volcanoes (Del Pezzo et al., 2006;De Siena, Thomas, & Aster, 2014;De Siena, Thomas, Waite, et al., 2014;De Siena et al., 2017;Prudencio, De Siena, et al., 2015;Prudencio, Ibáñez, et al., 2015) and tectonic contexts (Sketsiou et al., 2021).This paper presents novel 3D seismic attenuation and scattering models of Central-Southern Italy, interpreting them as images of the principal tectonic structures and fluid pathways.

Seismotectonic Setting
The Apenninic Chain is a Miocene-early Pliocene thrust-and-fold belt created by the convergence between Africa and Euroasiatic plates and the rifting of the Tyrrhenian Basin (Lavecchia, 1988;Patacca & Scandone, 2007;Patacca et al., 2008).The NE migration of the Chain is driven by W-dipping thrusts that superimpose inner stratigraphical domains on more external ones (Figure 1a).Late Pliocene-Quaternary high-angle SW-dipping normal faults dissect the pre-existing compressional structures (Barchi, 2010).These structures are responsible for devastating historical and instrumental earthquakes and represent upper-crust structures antithetic to significant systems of low-angle NE-dipping faults (Brozzetti, 2011;Brozzetti et al., 2009).Across peninsular Italy, the W-and E-dipping faults bound an upper-crust extensional seismogenic domain comprising most earthquakes of the Italian seismic catalog (Figure 1b, Lavecchia et al., 2021).The intra-Apennine extensional domain overlays a deeper compressive province developed at the hanging wall of the west-dipping outermost front of the Apennine compressional system (de Nardis et al., 2022;Ferrarini et al., 2021).Finally, the E-W striking transcurrent faults outcropping in the Gargano area and deepening toward the west (Argnani et al., 2009) define a strike-slip domain below the compressional one and contribute to the high seismic hazard of Central Italy.
Within the study area, the Quaternary extensional belt extends NNW-SSE across the Marche-northern Abruzzi sectors of the Central Apennines and propagates southwards in the WNW-ESE direction across the Molise-Northern Campania sector of the Southern Apennines (Figure 1a) (Carafa et al., 2020;Lavecchia et al., 2017Lavecchia et al., , 2022)).The Marche-northern Abruzzi sectors consist of many WSW-dipping high-angle normal faults, including the Fucino fault and the Paganica faults, responsible for Avezzano 1915 (Mw7, 10.000 casualties) and L'Aquila 2009 (Mw 6.3, 330 deaths) earthquakes, respectively.The outer front of the west-dipping system runs from the Gran-Sasso south-dipping structure to the WSW-dipping Morrone-Porrara alignment.The Marche-northern Abruzzi faults detach on an east-dipping basal detachment (Castaldo et al., 2018) referred to as Latium-Abruzzi extensional detachment, in Lavecchia et al. (2017) that deepens eastward to a depth of about 12-15 km.The Molise-Northern Campania extensional sector consists of eastward and westward fault systems developed in a prevailing WNW-ESE direction for a width perpendicular to a strike of about 40 km.The most external system is the eastward-dipping Isernia-Bojano fault alignment along the eastern side of the Matese Massif (Ferrarini et al., 2017).WSW-dipping normal faults dissect the western side of the Matese Massif.The Outer compressional Front has an arcuate eastward convex shape and is articulated in two minor arcs, the Abruzzo Citeriore Arc northward and the Frentani Arc southward (Ferrarini et al., 2021).The frontal thrust and its splay are buried under Plio-Pleistocene deposits, but they are evident in seismic lines and highlighted in morphotectonic studies.The Outermost Basal detachment and/or the Casoli and Bomba major inner splays are considered responsible for the most significant historical events of the area, which are the 1706 (Mw 6.8) and 1933 (Mw 5.9) earthquakes (DISS Working Group, 2015;Lavecchia et al., 2021).The westernmost sector of the study area comprises part of the Roman Comagmatic Province (Peccerillo, 1985) near the Tyrrhenian coast, more specifically, the area lacking volcanism between the Southernmost tip of the Albani Hills (Latium) and the Roccamonfina Volcano (Campania).The Province is the primary source of gas emissions in Central-Southern Italy (Chiodini et al., 2004;Frondini et al., 2019).

Data and Methods
Our starting data set comprises earthquake waveforms recorded within the study area (latitude: 41.1-42.5, longitude: 12.5-15.5)by the Italian seismic network in the period 2009-2023, corresponding to the stable configuration of the network.We selected available events having a magnitude range of 2.0 ≤ M L ≤ 4.5 and a depth  , 2007;Rovida et al., 2020).between 0 and 50 km.We inverted the coda-normalized direct-wave energies and mapped the peak delay times by using the MuRAT package (De Siena, Thomas, & Aster, 2014;De Siena, Thomas, Waite, et al., 2014;Reiss et al., 2022), which computationally implements coda-normalization and peak-delay techniques, providing a complete open-access tool for attenuation and scattering tomography.After downloading and checking the available waveforms, we used a final input data set that comprises 240 earthquakes and 7,287 waveforms (Figure S1 in Supporting Information S1).MuRAT also requires seismic velocities for tracing rays.However, no unique 3D velocity model exists for the area: we thus use as a reference the 1D model (Table S1 in Supporting Information S1) from Trionfera et al. (2019) extended over a grid spaced about 10 (horizontal) and 5 km (vertical).
Before applying MuRAT, we analyzed the frequency contents, the P-and S-arrival time distribution, and the start time of the coda in the SAC waveforms (the statistical analysis can be found in the Text S1 in Supporting Information S1).The pre-processing guided the selection of central frequencies (Figures S2 and S3 in Supporting Information S1) and other inputs for the attenuation and scattering tomography.We tested the stability of the attenuation and scattering measurements (Figure S4 in Supporting Information S1), explicitly checking that (a) the logarithm of peak delay increases linearly with source-station distance and (b) the logarithm of the total coda-normalized P-wave decreases with increasing travel time (Sketsiou et al., 2021).We used P-phases as a reference for the attenuation tomography because of the high number and quality of the detections for this phase; instead, we used S-waves for the peak-delay model due to the sensitivity to heterogeneities and the lack of tradeoff with the Vp/Vs structures (De Siena et al., 2016;Takahashi et al., 2009).
The total attenuation and scattering measurements have been performed at 3 and 6 Hz, respectively; they both have ray-dependent sensitivity.A standard regionalization procedure, which divides the volume into cells having the same values, is applied for the 3D peak-delay (PD) model (Gabrielli et al., 2023), while a different representation method is used for the Q model: MuRAT inverts coda-normalized direct-wave P-wave energies for the total attenuation (Q) (Del Pezzo et al., 2006;Sketsiou et al., 2021) and inverts the attenuation values along the seismic raypath.Checkerboard and spike tests for this method show differences in resolution across the entire region and for interpreted anomalies, respectively (Figures S5-S9 in Supporting Information S1).We realized an interpolated map of values for each depth layer using the Inverse Distance Weighted interpolation method.Only anomalies recovered by testing and that can be recognized in both the Q and PD models were interpreted in relation to existing geological features.

Attenuation and Peak-Delay Models
We provide new 3D models of total seismic attenuation (inverse total quality factor-Q model) and scattering (Peak-Delay models, hereinafter PD) for Central-Southern Italy (Figure 2).The horizontal slices of the attenuation model show four significant anomalies (labeled A-D in Figure 2).Their dimensions are comparable with the average resolution of the models (about 20 km), and the synthetic tests demonstrate they are resolvable at different frequencies (Figures S5-S12 in Supporting Information S1).
Along the Apennines Chain, at depths of 5.6 km, a nearly continuous high-attenuation zone extends NW-SE from L'Aquila to Benevento.This zone is identified in Figure 2b and exhibits two prominent anomalies labeled A and B. Anomaly A spans approximately 110 km (NNW-SSE) and shows a complex pattern (a1, a2, a3 in Figures 2b  and 2c).Feature a1 is in the southern Apennines and extends primarily along the Matese Massif and the Volturno Plain, bordering the volcanic area of Roccamonfina to the southwest.The area is partially characterized by high CO 2 emissions at the surface (Frondini et al., 2019;Vitale et al., 2023) (Figure 3).Anomalies a2 and a3 mainly extend in the southern Abruzzo region, along the Apennine axis (a2) and eastward, on the foothills (a3).Anomaly B, located in the Central Apennines near the Rocca di Mezzo fault system, extends for ∼20 km in the SW-NE direction (Figures 2a-2c).It is south of L'Aquila town, which experienced a long and destructive seismic sequence in 2009.The spatial distribution of the seismicity associated with this sequence abruptly interrupts North of the B anomaly (Figure 2b).The third well-resolved shallow high-attenuation volume (C) is located on the eastern side of the study area, between Benevento and Campobasso (Figures 2b and 2c).At a depth of approximately 12 km, anomaly D elongates NW-SE between Benevento and Isernia (Figures 2d and 2e).This anomaly shows attenuation values greater than 0.015 and can be divided into two separate areas: the southernmost volume north of Benevento (d1) and the one beneath the Matese Mountains (d2).The seismicity between 8 and 13 km is located in the low attenuation region between d1 and d2 (Figures 2d and 2e).
Anomalies A, B, C, and D are roughly identified in the PD model (Figures 2c and 2e), but their boundaries are better defined in the Q model.Their shape could differ due to differences in regionalization and grid-based inversion methods.Therefore, we used total attenuation to reference their boundaries and shapes.Also, in the PD model and at shallow depths, anomalies A and B show high scattering values interspersed with lower to negative values (Figure 2d).The external shallow (C) and 12-km-deep (D) anomalies appear as high-scattering zones (Figures 2c  and 2e).The Q and PD maps agree on the sharp variation between A and B anomalies.Low-to-high attenuation and negative-to-positive PD characterize the area, including the fault segments of Morrone-Porrara and Fucino.

Control of Tectonic Structures on Seismic Attenuation and Scattering
The most likely cause of high seismic attenuation and scattering volumes are rock fracturing and morphological variations (Gabrielli et al., 2023;Napolitano et al., 2020).Our results reveal significant attenuation and scattering values within both extensional and contractional domains (Figures 1a and 2).These attenuation maps exhibit a  (Calamita et al., 2009;Ferrarini et al., 2021;Patacca et al., 2008), which would intercept the high-attenuation volume at about 6 km (assuming a mean dip angle of 20°).Its WSW-ENE trend is replicated in the north by anomaly B, perpendicular to the primary NW-SE trend of the central Apennine fault systems, located between the west-dipping Paganica and Salto Valley-Velino fault systems.B comprises three minor WSW-dipping en-echelon faults bounding a small Quaternary basin; the envelope of these faults appears to control the anomaly.
The C anomaly and the southern N-S trending segment of the Frentani Arc of the Outer Thrust System follow similar trends (Figures 2a-2c).The Q image shows an arcuate trend following the thrust trace, but its shape is unclear in the PD model.Here, we interpret the anomaly as associated with the southern portion of the Frentani Arc.The peak delay maps strengthen and confirm the results obtained by inverting for the total attenuation.The two main high attenuation volumes at shallow depths (A, C) correspond to regions with dense and complex fault patterns able to release significant earthquakes.These fractured shear zones trap seismic waves, dramatically increasing peak delays while waves propagate within their fractures (King et al., 2022(King et al., , 2023)).
The spatial relation between total and scattering attenuation suggests a primary contribution of strain release to energy loss; even more relevant appears the link between low scattering and total attenuation and zones of stress accumulation (Gabrielli et al., 2023;King et al., 2022).While low scattering attenuation values are linked to compact, undamaged rocks, recent studies demonstrate their ability to define aseismic, locked volumes most likely to generate large earthquakes (Gabrielli et al., 2023), justifying the result as caused by a deformation of the heterogeneities due to the high-stress level in the proximity of the faults (King et al., 2022).A WSW-ESE-trending aseismic low 1/Q and low PD gap extending continuously for ∼120 km transversely to the Apennine fault system, from the Latium Volcanic Province to the Maiella Massif, separates anomalies B and a3 (Figures 2b and 2c).In the past, several destructive earthquakes struck the area, like the 1706 and 1933 events (Io greater than IX, Guidoboni et al., 2019).We propose that this low-attenuation and low-scattering barrier for earthquakes comprising exposed normal faults (Frepoli et al., 2017;Romano et al., 2013) is a locked zone where seismic stress accumulates, prone to release strong earthquakes.In particular, the Morrone fault system is known as a relevant seismic gap in the central Apennines (Chiarabba et al., 2009;Galadini & Galli, 2000;Gori et al., 2011;Romano et al., 2013;Vignaroli et al., 2022).
At a depth of about 12 km, the only anomaly with very high attenuation values is D, which is potentially controlled by the east-dipping structures bordering the Volturno plain.This anomaly is better related to fluid distribution and gas emissions across Southern Italy.

The Role of Fluids
Seismic absorption is better suited than scattering or total attenuation when imaging fluid reservoirs and molten bodies (Napolitano et al., 2020;Reiss et al., 2022;Sato et al., 2012).However, diffusive coda waves present absorption sensitivity kernels that depend heavily on source and station locations (Del Pezzo et al., 2018); differently from small volcanic and tectonic regions, where seismicity can diffuse over long periods, at our scale 3D coda-dependent absorption tomography provides much lower resolution than ray-dependent imaging.After using MuRAT to obtain absorption tomography models, we realized these were also unstable, that is dependent on few available measurements in parts of the models.Therefore, we grounded our interpretation on the spatial relation between total and scattering tomography anomalies, as fractures strongly contribute to fluid migration and storage.
Additional factors such as rock composition, mineralogy, and the presence of fluids can also play a role in influencing the attenuation of seismic waves.The westernmost sector of the study area comprises the Southern tip of the Roman Magmatic Province.This sector is characterized by diffuse volcanic phenomena such as CO 2 degassing (Frondini et al., 2019;Vitale et al., 2023) (Figure 3).Here, a high-attenuation/high-scattering volume coincides with the Albani Hills area, underlying the well-known correlation between seismic attenuation and fluids sources in volcanic systems (Grab et al., 2017;Hudson et al., 2023;Schurr et al., 2003).Its features are unresolved (Figure S9 in Supporting Information S1), but values are similar to those of B (Figure 2), which comprises a Quaternary basin.While portions of B present artifacts due to the source-receiver distribution (Figure S9 in Supporting Information S1), it also corresponds to the volume interested by the 2009 sequence: the fluids that permeate fractures likely increase attenuation.We propose this transversal fracturing zone as a reservoir for the fluids involved in the seismic sequence (Martinelli et al., 2020;Quattrocchi et al., 2011;Terakawa et al., 2010).Nevertheless, our models have insufficient resolution to supply the exact shape of this reservoir.
A similar relation between fluids and seismic attenuation and scattering is proposed in the Southern Apennines, where D (d1, d2) dominates between 10 and 15 km (Figures 2d and 2e).In the PD map, this feature propagates mainly SW.It reaches higher scattering values in correspondence with the Campanian Volcanic Province and the associated large CO 2 flux (Figure 3d).The high 1/Q instead borders the high CO 2 flux zone (Figure 3c).The position of the anomaly relative to gas fluxes suggests a deep reservoir separated from the volcanic area of the Campanian Volcanic Province.The Volturno Plain normal fault alignment (Figure 2a) could represent the barrier that separates the two fluid sources, avoiding their mixing.The reservoirs created by the Campanian Volcanic Province degas through tectonic and volcanic structures, while d1 and d2 appear trapped at depths of about 12 km.These fluids can partially fill the tectonic fractures following the principal W-and E-dipping normal faults; consequently, they can reach shallower depths and contribute to the high-attenuation, high-scattering values of the wider A anomaly.This prominent feature is the first evidence of a fluid-saturated fracture zone that likely sources fluids to the fracturing system crossing the Southern part of our study area, and that could have a crucial role in the seismotectonic assessment and the earthquake triggering.Our results are consistent with the analysis of Di Luccio et al. ( 2018), which provided additional insights into the genesis of seismic events in this 10.1029/2023GL106074 8 of 11 area, attributing them to tectonic loading, over-pressurized CO 2 reservoirs at upper crust depths, or intrusive episodes of magmatic sources at greater depths.A recent paper by Vitale et al. (2023) also confirmed the presence of deep fluids in the area, unveiling four gas vents around the Matese Massif area (Figure 3).

Conclusions
The lack of a homogeneous geophysical model for Central-Southern Italy significantly limits the understanding of shallow and deep crustal structures and their hazard assessment.In this paper, we address filling this gap by obtaining 3D total attenuation and scattering models for the area.The two independent parameters reveal similar wide-scale features, reinforcing their reliability and allowing a multiparametric interpretation.Major tectonic structures control the distribution of high-attenuation volumes.In particular, high attenuation in Southern Italy is related to the W-dipping Boiano, Monte Greco-Barrea, and the E-dipping Volturno-Benevento Fault Systems.We also detect the Southern prosecution of the Frentani Arc as an arcuate high-attenuation zone.
High scattering and attenuation characterize a deeper fluid-saturated volume separated from the Campanian Magmatic Province by the Piedimonte Matese fault system.We propose this reservoir as the source region for shallower fluid circulation through the tectonic fracture network that interests the South of the inspected area.In the Central Apennines, diffuse low values of PD and 1/Q mark seismic stress accumulation in the area of the Morrone-Porrara fault alignment, where they coincide with a seismic gap.At the southern and northern borders of this volume, high-attenuation anomalies rotate WSW-ENE perpendicular to the Apennine axis.These signatures lead us to propose transversal fluid-filled volumes that "stop" earthquake propagation toward the North and South by reducing the normal stress on the faults and favoring an almost continuous stress release.Due to the stress accumulated, large earthquakes could strike the aseismic low-attenuation and scattering gap in the future.

Figure 2 .
Figure 2. (a) Summary of the main tectonic features in the area; faults colors are the same as in Figure 1.(b-e) Horizontal maps of total attenuation parameter (Q) (b, d) and scattering proxy (peak-delay [PD])(c, e).Q parameter refers to P-waves, and PD refers to S-waves.The rows show results at depths of 5.6 and 11.7 km, respectively.The seismicity distribution of Figure1is projected on each map with a semi-width of 2.5 km: earthquake ranges are 3-8 km (5.6 km, (b, c)) and 8-13 km (11.7 km, (d, e)).Black dashed lines and labels are the anomalies described in the main text.

Figure 3 .
Figure 3. Map of surface CO 2 flux in Central-Southern Italy interpolated from the degassing measurement in correspondence with springs (Frondini et al., 2019).White dashed and dotted lines are the contours of anomalies in the Q(p) and peak-delay (PD) models.The CO 2 flux (Frondini et al., 2019; Vitale et al., 2023) is compared with the anomalies from the Q and PD models at 5.6 km (a, b) and 11.7 km (c, d).