3-D attenuation image of fluid storage and tectonic interactions across the Pollino fault network

P. Sketsiou ,1 L. De Siena ,1,2 S. Gabrielli 1,3 and F. Napolitano4 1School of Geosciences, University of Aberdeen, Aberdeen AB24 3UE, Scotland, UK. E-mail: panayiotasketsiou@gmail.com 2Institute of Geosciences, Johannes Gutenberg University, D-55128 Mainz, Germany 3Instituto Nazionale di Geofisica e Vulcanologia (INGV), 00143 Rome, Italy 4Dipartimento di Fisica “E.R.Caianiello”, Universitá degli Studi di Salerno, 84084 Fisciano (SA), Italy


I N T RO D U C T I O N
Seismic waves lose energy during propagation in heterogeneous Earth media. They lose amplitude in space and time due to seismic attenuation, whose description is central when modelling seismic wave propagation both at the laboratory and field scales (Müller et al. 2010;Cormier 2011;Sato et al. 2012). The total quality factor (Q) measures the anelastic attenuation of coherent waves. It is defined as the fractional energy lost per cycle and controls the decay of the energy density spectrum with lapse time from the origin time of the earthquake (Cormier 2011). While in many seismological studies the physical processes measured by Q (e.g. scattering and absorption) are considered a hindrance and removed, the attenuation of body, surface and coda waves has become a 536 C The Author(s) 2021. Published by Oxford University Press on behalf of The Royal Astronomical Society. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
The most widely used method to estimate the attenuation of seismic waves is the so-called t * method (Romanowicz 1998;Schurr et al. 2003), where t * represents the cumulative attenuation along a ray path (Kanamori 1967). The t * method measures the high-frequency spectral decay rate above the corner frequency (dependent on earthquake magnitude) for each waveform: e −ωt * /2 (Cormier 1982;Bindi et al. 2006), where ω is the angular frequency. While its application is easy and justifiable in the presence of gradually varying heterogeneity, the assumptions underlying the method are often unfulfilled if there are strong velocity variations, interference between coherent waves, and complex source and site spectra (Del Pezzo 2008;De Siena et al. 2009;Parolai et al. 2015). The method does not allow a frequencydependent reconstruction of structures, even if these certainly interact differently with seismic waves depending on sampling frequency. Aki & Chouet (1975) and Rautian & Khalturin (1978) developed the coda normalization as an alternative method to measure Q of S waves (Q S ) within highly heterogeneous lithospheric media. Yoshimoto et al. (1993) extended the method to measure the attenuation of P waves (Q P ), normalizing P-wave energy by S-coda energy. The effect of different radiation patterns on coda-normalized measurements is negligible if the data set consists of events with a wide range of epicentral distances that allow to average over many focal mechanisms (Yoshimoto et al. 1993). Even considering singlestation measurements, radiation patterns only have secondary effects if the window taken to measure direct-wave amplitudes is long enough. Takemura et al. (2009) studied the distortion of the S-wave radiation pattern of small-magnitude earthquakes due to small-scale crustal heterogeneity. They concluded that the S-wave radiation distributes quasi-isotropically for high-frequency waves after a few mean-free times from the origin-time of the earthquake. Obviously, when enlarging the window used to measure P-wave amplitudes, the risk is to include spurious surface phases, finally losing focus on the infinite-frequency seismic ray. A preliminary measurement of spatially varying coda properties is also necessary (De Siena et al. 2014a).
The coda normalization method is a standard for measuring and mapping the average attenuation properties of the lithosphere. Examples comprise the comparison of shear wave attenuation in different areas of the crust (Frankel et al. 1990) and collision zones (Sharma et al. 2007;Parvez et al. 2012). Del Pezzo et al. (2006 first developed an attenuation tomography method based on coda normalization (CN) to image Mount Vesuvius. They measured Q S along rays traced in a velocity model to invert for the 3-D attenuation structure of the volcano. De Siena et al. (2009) benchmarked the method against t * , obtaining a multiple resolution 3-D model of both Q S and Q P at the same volcano, using a direct-wave time window of 3 s to remove radiation pattern effects. At Campi Flegrei caldera, where topographic effects are marginal, the window length was reduced to 2 s (De Siena et al. 2010). Since then, body-wave attenuation models based on the CN method have been applied at crustal scale across the world, becoming a standard in volcanoes (Matsumoto et al. 2009;De Siena et al. 2014b;Prudencio et al. 2015a, b;De Siena et al. 2017;Nazemi et al. 2017;Prudencio & Manga 2019).
While the coda-normalization method has been used to study the average attenuation properties of fault zones (Akinci & Eyidogan 1996) a 3-D tomographic inversion has never been applied in this setting. De Siena et al. (2014a) have published MuRAT, a MATLAB R code to obtain frequency-dependent Q P and Q S measurements using the CN method and invert for the P-and S-wave attenuation structures of highly heterogeneous regions. In its 2.0 version, the code allows to estimate coda attenuation before normalization, reducing the unknowns in the inversion.

S T RU C T U R E A N D DY N A M I C S O F T H E P O L L I N O FAU LT N E T W O R K
The data in the present work have been recently used to map the 2-D frequency-dependent spatial variations in scattering and absorption across the Pollino fault network (Napolitano et al. 2020a). This fault network comprises low-angle, E-and NNE-dipping faults, as well as antithetic, high-angle, SW-to WSW-dipping faults (Brozzetti et al. 2017). The area is historically recognized as a seismic gap in the Apennines chain (Passarelli et al. 2015). Regardless, there is a record of historical earthquakes with magnitudes between 5.2-6.0 dating back to the seventeenth century (Tertulliani & Cucci 2014). The Pollino seismic sequence has occurred between 2010 and 2014 within the gap and with a highly variable rate.
The activity started in Fall 2010, with a small cluster of events, which faded out in Spring 2011 and started again in Fall 2011, located in the SE part of the Mercure Basin. In May 2012, an event with M w = 4.3 occurred to the east of the first sequence, with a classical main shock-aftershock sequence (hereafter "cluster 2", as labelled in Fig. 1). In July 2012, the activity restarted in the first location, in the SE of the Mercure Basin, lasting several months, with a peak M w = 5.1 normal faulting event occurring on 2012 October 25. Although the seismicity gradually faded by Spring 2012, the seismicity rate was higher than usual until 2014. The seismic activity in this location will be hereafter called 'cluster 1'. These sequences comprised thousands of low-magnitude earthquakes but no main shock-aftershock succession (Passarelli et al. 2015;Totaro et al. 2015;Ferranti et al. 2017). Such features, together with the occurrence of slow-slip earthquakes (Cheloni et al. 2017), have supported the identification of Pollino as a region of slow deformation. This slow-slip event, equivalent to an M W ∼ 5.5, has been detected between August 2012 and March 2013, and located in the same seismogenic volume of the first cluster. Increased pore pressure and, more generally, pressurized fluids modulating stress along the fault, are the most likely sources of slow-slip events and aseismic transients seismicity (Passarelli et al. 2015;Cheloni et al. 2017).
A frequency-dependent scattering and absorption mapping has provided the spatial measurements of coda attenuation necessary for the application of MuRAT2.0 (Napolitano et al. 2020a). This work has brought forward the hypothesis that lateral fluid migrations from the Castrovillari basin, SE of the fault network, is the most likely source of historical seismic events (Napolitano et al. 2020a;Sketsiou et al. 2020). Their high-scattering and high-absorption patterns shift towards NW following the strike of NW-to-SE-trending faults. The SE-NW migration stops at the boundary of the eastern unit of the Apennine Platform (Brozzetti et al. 2017) at intermediate frequencies. It marks the faults where the Pollino seismic sequence nucleated in the highest frequency range. Despite the increasing amount of shallow geological and geophysical information (Passarelli et al. 2015;Totaro et al. 2015;Brozzetti et al. 2017), leading to new models and hypotheses of the deep fault structures (Ferranti et al. 2017), deep 3-D geophysical models of the area are still lacking.

DATA
The data set consists of 102 earthquakes (for a total of 702 waveforms) that occurred between 2010 and 2014 in the seismogenic volume responsible for the Pollino seismic swarm-like sequence (Passarelli et al. 2015;Napolitano et al. 2020b) and nearby, as shown in Fig. 1. The waveforms have been recorded at 18 seismic stations operated by Istituto Nazionale di Geofisica e Vulcanologia (INGV), Universitá della Calabria (Margheriti et al. 2013) and GFZ, Potsdam, Germany (Passarelli et al. 2012). The strongest events of this data set (M L > 2.7) have been used to study the source mechanism and the development of the swarm (Totaro et al. 2015) and to constrain the geometry of the faults inferred by Brozzetti et al. (2017). Most of these waveforms have been used to assess ground-motion amplifications (Napolitano et al. 2018) and separation of scattering and absorption contribution to the total attenuation of coda waves (Napolitano et al. 2020a;Sketsiou et al. 2020). To study the P-wave attenuation, we selected the vertical component of the Pollino waveforms, characterized by source-receiver distances between 1 and 40 km, by local magnitude between 1.8 and 4.3, and depths between 1 and 32 km. The equipartition between P and S waves (Napolitano et al. 2020a) allows us to use the vertical component to compute coda waves at late lapse times, where the diffusion regime is reached (Sketsiou et al. 2020). The selected area, [lon: 15.4-16.4 (111 km); lat: 39.6-40.4 (70 km); depth: 20 km] was divided in blocks of 0.1 • , a suitable size to solve a 3-D problem. For ray tracing, we interpolated the velocity model by Barberi et al. (2004) between depths of 0 and 70 km, tracing rays with a ray-bending approach (De Siena et al. 2014b).
We performed two different tomographic analyses, corresponding to two data sets. Data set 1 comprises all the waveforms described above to invert for the entire volume under consideration. Data set 2 is a selection of 647 waveforms chosen to provide a higher-resolution image of the zone where the Pollino swarm developed: [lon: 15.8-16.2 (44 km); lat: 39.8-40.05 (21 km); depth: 12 km].

Assumption of frequency dependence
The (in)dependence of Q on frequency has been a debated question for decades and is a central assumption for any imaging technique using attenuation as attribute (i.e. attenuation tomography). When Q is considered independent of frequency (Knopoff 1964), a valid assumption for measuring attenuation, for example, in rock-physics (Tisato & Madonna 2012), the frequency-shift method becomes the ideal technique to image Earth media (Quan & Harris 1997). Still, the complex heterogeneous attenuation of the lithosphere causes simultaneous phenomena which span, among others: (1) the enhancement of high-frequency near-receiver surface waves due to scattering (Campman et al. 2005;Gabrielli et al. 2020); (2) anisotropic attenuation due to fractures (Chapman 2003); (3) the wave-induced flow acting at frequencies lower than 1 Hz (Pride et al. 2004;Müller et al. 2010). Most researchers agree that Q measurements are dependent, if not proportional, to frequency, at least at field scale and above 1 Hz (Fedotov & Boldyrev 1969;Rautian & Khalturin 1978;Aki 1980;Chapman 2003;Li & Lu 2010). The following imaging is thus performed at multiple frequencies. In zones of high resolution, frequency-dependent images can be interpreted in terms of scale-variations in geological structures and tectonic activity (Sato et al. 2012).

The coda normalization equation
Seismic wave amplitudes can be modelled as the convolution of the source, path, site, and instrument functions in the time domain (Del Pezzo et al. 2006). The spectral amplitude of direct S waves : where f is the frequency, r is the source-receiver distance, t c is the coda-wave central lapse time, R θφ is the source radiation pattern (θ and φ the azimuth and take-off angle for a source-receiver ray), S S (f) is the source function, γ the geometrical spreading exponent, I(f) the known instrumental response, G(f, ) is the site amplification factor (with being the incident angle of the ray at the station), Q S is the S-wave quality factor and v S the average S-wave velocity in the medium (Yoshimoto et al. 1993).
At late lapse times, where the diffusion regime is reached, the spectral amplitude of coda waves can be written as ( The coda amplitude must be independent of the hypocentral distance (Calvet & Margerin 2013): for Pollino, this was proven by previous works (Napolitano et al. 2020a;Sketsiou et al. 2020). The attenuation of coda waves can then be expressed as (Aki & Chouet 1975): In this equation, n is the envelope spectral decay, Q c is the coda quality factor and A c (f, t c ) does not include the effect of the source radiation pattern. Eq.
(1) can be divided by Eq.
(2) to normalize the energy of the spectral amplitude of direct S waves using the spectral amplitude of coda waves: where the source function S S (f) and the instrumental response I(f) contributions disappear. If the length of the direct-wave window is chosen appropriately, the contribution of the source radiation pattern R θφ is negligible (De Siena et al. 2009). Since early coda consists of randomly scattered waves, the coherency and the source radiation pattern is eventually lost (Takemura et al. 2009). De Siena et al.
(2010) demonstrate that a coda window of 2 s is sufficient to make the radiation pattern quasi-isotropic in a volcanic caldera. Still, the window length must be carefully chosen at each frequency to avoid near-receiver onset of surface waves (Gabrielli et al. 2020). At fixed frequency band and starting lapse-time for coda windows, P(f, t c ) is assumed to be constant in the standard codanormalization method (Del Pezzo et al. 2006;Sato et al. 2012). Taking the logarithm, eq. (4) becomes By pre-computing G(f, ) over many earthquakes we remove its dependency on the incident angle . Finally, dividing by π f, the coda-normalization equation is defined as where Yoshimoto et al. (1993) demonstrated that earthquakes within a small magnitude range show the same spectral shape for P and S waves. Hence, we can write and that Eq. (6) can be used for coda normalization using P-wave velocities (v P ): Both Eqs (6) and (8) are forward models that can be used in tomography to solve for quality factor variations (Del Pezzo et al. 2006;De Siena et al. 2010, 2014b. In the following, these forward models are used by MuRAT (Del Pezzo et al. 2006;De Siena et al. 2014a) in order to map Q at Pollino. The new version, currently integrated in MuRAT2.0 (https:// github.com/LucaDeSiena/MuRAT), pre-computes spatially varying coda attenuation (Sketsiou et al. 2020) and uses it to remove one of the unknowns (K c (f)) from Eqs (6) and (8).

Tomographic procedure
We follow the methodology developed by Del Pezzo et al. (2006) and implemented computationally by De Siena et al. (2014a). This methodology uses Eqs (6) and (8) to obtain a tomographic forward model. From the linear regression of Eqs (6) and (8), the average inverse total quality factor is Q −1 P,S (where the subscripts refer to either P-or S-wave amplitudes), the geometrical spreading (γ ) and K c (which is assumed to be constant for a fixed frequency band and lapse time) are obtained. As explained before, the only unknown in K c is the average coda quality factor Q c (Eq. 3), which is measured experimentally for each source-receiver pair.
This pre-calculation allows to define a new data vector (d k ) for each source-station path (k): The area under study is divided into a grid of M equally spaced nodes and the forward model is derived from Eqs (6) and (8) where k refers to the source-station path, B indicates the Bth of the M blocks that the kth ray crosses, s B is the slowness of segment of length l k B crossing the Bth block, and δ(Q B P,S ) −1 are the variations from Q −1 P,S in each block. Eq. (10) is solved using a zero-order Tikhonov regularization. The damping parameter can be chosen from the corresponding L-curve (Aster et al. 2013). The total quality factor in each block is obtained as

Average direct and coda attenuation values and dependence on frequency
We start with reporting the results obtained using Data set 1 (see Sketsiou et al. (2020) for details). The dependence of ln Q −1 P on the 540 P. Sketsiou et al.  Table 1. Best set of input parameters. t c and L W are the starting lapse-time and length of the coda window, respectively and L P is the length of the window for the body wave analysis. The window for the analysis of the body waves starts with the P-wave onset, and is chosen taking into account the t S − t P of the waveforms with the shortest epicentral distance. logarithm of frequency is shown in Fig. 2(a). The points were fitted by taking the logarithm of an equation of the form Q −1 Petrosino et al. 2008). Several parameters affecting the results such as the coda window onset t c and length L w , length of the P-wave window L P and the envelope spectral decay n have been tested (Havskov et al. 2016). The best set of these parameters are shown in Table 1. The values of Q P and γ obtained using Eq. (9) for frequencies between 3-18 Hz are shown in Table 2. The resulting equation is Q −1 P = (0.0011 ± 0.0008) f (0.36±0.32) . The p-value of the regression model, which tests whether this model fits significantly better than a degenerate model consisting of only a constant term, is equal to 0.343 which is higher than 0.05, meaning that the aforementioned linear model does not adequately represent the data. Additionally, we modelled the dependence of Q −1 c on frequency. The points were fitted with a power law as in the case of Q −1 P , resulting in the relationship Q −1 c = (0.0057 ± 0.0006) f (−0.55±0.05) . Table 2. Q −1 P and geometrical spreading γ measurements with their associated errors σ Q P and σ γ (one standard deviation) obtained from the linear regression at different frequency bands using the best set of parameters, shown in Table 1.
Data set 2 (Section 3), we calculated the same dependence of Q −1 P on frequency. This data set spans a depth range of just 12 km, half of Data set 1, and are constrained around the source region. With the same set of input parameters as described for Data set 1 (Table 1), the resulting relationship is showing the expected higher attenuation and frequency dependence. The p-value of the regression model is equal to 1.61 × 10 −6 , which is significantly lower than 0.05, meaning that the linear model represents the data well. However, the dependence of Q −1 c on frequency results in the relationship Q −1 c = (0.0057 ± 0.0006) f (−0.55±0.05) , which is identical to Data set 1. This is a strong indication that coda waves behave homogeneously in space regardless of the size of our study area at the chosen lapse times. The Q −1 0 and η values obtained are characteristic of stable continental crust (Data set 1) and high-seismicity zones (Data set 2), respectively (Sato et al. 2012).

Frequency-dependent tomographic results
We obtain several 3-D models of the Pollino fault network in different frequency bands and interpret them in terms of both resolved features and scale of the heterogeneity. We consider horizontal slices taken at 2.5 km, 7.5 km and 12.5 km in the three frequency bands (Figs 3(a)-(i)) as well as S-N, W-E and NW-SE vertical profiles crossing the area of the Pollino swarm (Supporting Information Fig. S1).
Given the velocities we used (Barberi et al. 2004) and the depths sampled by the earthquakes, we expect to resolve structures of the order of 10 km 3 across the medium at 3 Hz (Fig. 4, top row). The only resolved anomaly at a depth of 2.5 km at all frequencies is a low-attenuation pattern located just above the volume where the last Pollino swarm occurred (Figs 3(a)-(c)). At 7.5 km, most of the volumes west of 16.3 • longitude and centred at 7.5 km are resolved (Fig. 4, middle row). Here, high attenuation anomalies trending NW-SE (Fig. 3(d)) appear across the northern fault network, following historical seismicity (Fig. 1). The area of the Pollino swarm, centred at this depth around [16.0 lon, 39.9 lat], has average attenuation characteristics. At a depth of 12.5 km, the resolved region is between 15.9-16.3 lon and 39.8-40.04 lat (Figs 3(g), S1). Two high attenuation anomalies are located NW and SE of the shallower epicentral area of the Pollino swarm ( Fig. 3(g)), respectively. The NW high-attenuation anomaly is at the borderline of our resolving power (Fig. 4, bottom row) and possibly a ghost. The resolved SE high-attenuation anomaly corresponds to the oldest historical seismicity recorded in the area (Tertulliani & Cucci 2014).
A NW-SE vertical cross-section (Supporting Information Fig. S1) shows that the low-attenuation volume topping the epicentral area extends from the surface down to a depth of 5 km above the Pollino swarm. The N-S and W-E profiles, taken at 16.01 lon  . Results of the checkerboard tests for the three depth slices at 2.5 (top row panels), 7.5 (middle row panels) and 12.5 km (bottom row) depths, at three central frequencies: 3 Hz (second column), 6 Hz (third column) and 15 Hz (last column). The first column shows the checkerboard inputs. 542 P. Sketsiou et al. and 39.92 lat, respectively (Supporting Information Fig. S1) show that this medium-to-low attenuation volume characterizes the whole southern portion of the study area. While the rupture zone for the M5 Pollino earthquake is sandwiched at the top of the two deeper high-attenuation anomalies in the NW and SE, only the SE anomaly (1) is resolved and (2) is related to historical seismicity. The only resolved high-attenuation anomaly at the average depth of the swarm (7.5 km) is north of it (Figs 3(d)-(f) and Supporting Information Figs S1(d)-(f)).
The 6 Hz model confirms the broadest features visible in the 3 Hz maps (Figs 3(b),(e),(h)): (i) a wide low-attenuation volume extending above and to the south of the Pollino hypocentral area down to 5 km; (ii) the NW-SE trend of high attenuation volumes centred at a depth of 7.5 km and following historical seismicity; (iii) two high-attenuation volumes sitting just below the hypocentral area (12.5 km depth), although again the northern anomalies are unresolved; (iv) the high-attenuation volume at 7.5 km depth, located just north of the hypocentral area showing average attenuation characteristics.
The 15 Hz maps provide a more resolved reconstruction of all the high-and low-attenuation anomalies listed above (Figs 3(c),(f),(i)). However, the high-attenuation anomaly north of the hypocentral area at 7.5 km now extends southwards along the NS-trending faults west of the swarm area.
The resolution tests at 2.5 km depth show that the only resolved area is directly above the swarm for Data set 1 (Fig. 4, top row). At 7.5 km, the majority of our study area is resolved, up to 16.3 • lon using Data set 1 (Fig. 4, middle row). Using Data set 2, the majority of the area between 15.88 • and 16.14 • lon is resolved at this depth (Fig. 5). At 12.5 km depth, the resolved area for Data set 1 is between 15.9-16.3 • lon and 39.8-40.04 • lat (Fig. 4, bottom  row).
For Data set 2 the resolution is limited both above and below 7.5 km due to the lower ray-crossing (Figs 5(b),(e),(h)). Regardless, they provide more details in the area of the seismic sequence than Data set 1 around a depth of 7.5 km (Figs 5(b),(e)). The highattenuation anomaly previously comprising the Mercure basin now extends across the western side of the cluster 1 events and stops in the NW and W due to a sharp low-attenuation interface, at the boundary of our lateral resolving power (Figs 5(c),(f),(i)). The lowto-average attenuation anomaly previously characterizing the entire seismic sequence now comprises only the eastern part of the cluster 1 events, as well as the aseismic zone between clusters 1 and 2. Fig. 6 is a simplified sketch of our tomographic model and includes geophysical and geological results used in the following interpretation. We primarily refer to this sketch while mentioning the tomographic sections when discussing resolution and reliability.

Frequency dependence of seismic attenuation measurements
Using Data set 1, there appears to be no dependence of Q −1 P on frequency ( Fig. 2 (a)). Considering the error on each value, we see that the trend is approximately constant, showing nearly no dependency on frequency. The reduction of the study area with Data set 2 led to focusing on the fault networks and, as a result, the Q −1 P dependence on frequency is significant. The Q 0 and η values obtained are indicative of tectonically active areas (Padhy 2005;Sato et al. 2012). The Q −1 c0 and η values obtained using Data sets 1 and 2 are the same: they indicate that coda waves behave similarly regardless of the size of our study area and they are characteristic of continental areas with high seismicity (Sato et al. 2012).

Low-permeability shallow volumes topping lateral seismic migration
The low-attenuation volume topping the hypocentral area of the cluster 1 events (depth of 2.5 km) corresponds to a relative absence of extended faulting and, consequently, lack of historical and more recent seismicity (Fig. 6(a); Tertulliani & Cucci 2014;Ferranti et al. 2017). This volume seems able to stop the vertical rupture of the western fault during the Pollino sequence. As fluids are the likely cause of swarms and aseismic transients (Passarelli et al. 2015;Cheloni et al. 2017;Napolitano et al. 2020a), the volume could mark semi-impermeable formations. The presence of such an impermeable layer at 4.5-5 km depth and lat >39.92 • is supported by the highly detailed relative location of clusters having similar waveforms, obtained by Napolitano et al. (2020b) and shown in its figs. 10 and 11. The hypocentres migrated horizontally (between 4.5 and 6 km depth) from the middle of the fault network towards North between Nov. 2011 and Feb. 2012, just below the interpreted impermeable formation. The re-located hypocentres started propagating vertically between Sept. 2012 and June 2013, only after the seismicity migrated towards the South over the boundary of the low-attenuation anomaly (lat <39.92 • ). The results of these two independent analyses obtained using different data sets, consistently point to a mechanism of lateral fluid migration and, possibly, excess of pore fluids (Passarelli et al. 2015), as the cause of the Pollino seismic sequence. The distribution of deep-sourced springs shows that these deep-sourced fluids find their way to surface along the SW boundary of the formations (Gaglioti et al. 2019), in spatial relation with the main aseismic transient, synchronous to the swarm (Cheloni et al. 2017) (Fig. 6).
At these depths, impermeable formations would correspond mostly to the transition from the shallower ductile-behaving Apennines imbricate wedge and the brittle Apulia sedimentary crust (Ferranti et al. 2017). According to Brozzetti et al. (2017), the location comprises a part of the Mercure basin, a part of a Liguride unit and a part of the Apennine platform western unit. The first ∼ 2 km layer of sediments in the area (limestones and dolomites) have thus likely high permeability.

Deep fluid reservoirs and low-attenuation, low-seismicity gap between the two clusters
All the higher resolution maps (Fig. 5) show that a wide highattenuation volume develops at a depth of 7.5 km across the Mercure basin, and comprises the western and central part of the cluster 1 events (Figs 5 and 6). Constrained in the north, west (Lauria) and east (Pollino Range) by the low-attenuation eastern unit of the Apennine Platform (Knott 1987;Giaccio et al. 2014;Brozzetti et al. 2017; Fig. 6), the high-attenuation volumes are likely representing the seismic signature of a deep fluid reservoir in the area ('Mercure reservoir'). These deep fluids find their way to the surface across the boundary between the reservoir and the Lauria mountains, the primary region of deep-sourced springs in the study area (Margiotta et al. 2012;Canora et al. 2019).
A low-attenuation volume develops inside the low-seismicity gap between cluster 1 and cluster 2 events (Figs 5 and 6(b)). It separates the Mercure reservoir from a second high-attenuation anomaly, east of the cluster 2 events and underlying the Castrovillari basin. This high attenuation volume dips towards SE, below 10 km ( Fig. 6(b)). The low-attenuation gap underlies (likely young) smallscale cross-faulting (Brozzetti et al. 2017) and comprises patches of high attenuation that increase in size at higher frequencies. There is no clear indication of consistent high-attenuation pathways between the Mercure and Castrovillari reservoirs, and these patches are below our resolution power. However, the attenuation values in this volume are higher than those in the surrounding regions (parts of the eastern units of the Appenine platform) (Fig. 6). These observations are hardly explained by the presence of a barrier for fluid migration similar to that inferred under the Lauria Mountains (Fig. 6). They better agree with the hypothesis that the volume comprises a system of asperities (Passarelli et al. 2015) that can break in the low-attenuation matrix.
The attenuation structure obtained using Data set 2 was compared with the velocity model of Totaro et al. (2014), which is different from the 1-D velocity model of Barberi et al. (2004) we used for our tomographic procedure. In their paper, Totaro et al. (2014) studied the crustal V P , V S and V P /V S structure of the southern Apennine-Calabrian arc border between 3 km and 30 km. For our study, the original model was interpolated between 3 km and 12.5 km and cross-sections at 3 km, 7.5 km and 12.5 km were plotted for our area of interest (Figs 5 and S3). In Fig. 5 positive and negative variations of attenuation are spatially correlated with positive and negative variations of V P /V S . The low V P /V S and attenuation support the interpretation of the low-seismicity gap as a poorly permeable system that slows down fluid migration.
The impermeable formations capping the cluster 1 events area and this deep low-attenuation volume depict the ideal setting for the development of transient slip events due to fluid redistribution and migration within the Mercure reservoir and adjacent faults (Passarelli et al. 2015). The aseismic transient identified during the seismic sequence by Cheloni et al. (2017) is located just under the impermeable formations. The transient develops NW across the reservoir, consistent with the transtensional stress field acting in the southern part of the Mercure Basin, associated with NNW-SSE normal faulting (Passarelli et al. 2015).

Fluid sources and earthquake-inducing fluid migration
Our model depicts a closed reservoir system lacking highattenuation fluid sources that could create pressure directly below 544 P. Sketsiou et al. Figure 6. Schematic interpretation of the attenuation model, as described in the text, illustrating the main findings. In panel (a), we show a horizontal section of the region with the small grey dots representing the seismicity, and including the location of the 1998 Mercure earthquake (Brozzetti et al. 2009). The dashed orange shapes represent high attenuation regions and the dashed blue shapes represent shallow low-attenuation volumes. The solid blue lines represent low attenuation regions related to mountain platforms. Three areas of increased stress before the swarm are marked with a white cross as well as two areas of aseismic transient. Surface springs with deep sources are marked with green circles. Panel (b) is an E-W cross section at 39.9 • latitude. The seismicity of the region at depth, the low-and high-attenuation volumes and the aseismic transient are shown, as described before for panel (a). We also show the deep low-seismicity and low-attenuation volume in light purple. The dashed black line marks the borderline of our resolution and the solid black line marks the limit of the tomographic model. the cluster 1 events (Fig. 3). Our models lose resolution at depths below 10 km, especially when using Data set 2 (Fig. 5) and are clearly more affected by the interpolation of the broader velocity model (Barberi et al. 2004), necessary to obtain high-resolution attenuation maps. We observe deep anomalies (12.5 km) west of the seismic events, corresponding to old extended NS-trending faults (Brozzetti et al. 2017;Figs 3(g)-(i)). These anomalies are unresolved and, given their aseismicity (both recent and historical), they are likely ghosts produced, for example, by the shallower Mercure reservoir. Despite the uncertainties, at least at the wider scale (Data set 1- Fig. 3), we are able to resolve a single deep high-attenuation volume that shows the extension to depth of the Castrovillari reservoir, SE of the fault networks (Figs 3 and 6).
The absence of high-attenuation volumes sourcing fluids just below the cluster 1 events (Figs 3(d)-(i) and 5) is expected given the structural characteristics of the area (Brozzetti et al. 2017;Ferranti et al. 2017) and the evidence of historical (Napolitano et al. 2020a) and recent (Passarelli et al. 2015;Napolitano et al. 2020b) lateral fluid migration. Historical earthquakes correlate with a general shift of fluid-related, high-absorption anomalies from the Castrovillari basin towards the NW, stopping at the foot of the Lauria mountains (Napolitano et al. 2020a). Even today, the deep volumes below the Castrovillari basin seem the most likely source of deep fluids. Here, at depths between 9 km and 14 km, Ferranti et al. (2017) speculate the existence of the ductile Verrucano (s. l.) layer, a siliciclastic and lower-velocity layer that rises above its average depth in the eastern part of our model. Similar ductile characteristics are detected as high-attenuation anomalies when mapping volcanoes with seismic attenuation (De Siena et al. 2017), thus they could affect imaging of the deeper volumes.
The trace of the Mercure reservoir at 7.5 km depth is close to the hypocentral location of the Mercure basin seismic sequence in 1998 (Brozzetti et al. 2009) (see main shock in Fig. 6). Our tests ( Fig. 4 and Supporting Information Fig. S2) show that the interface under the Lauria mountain is a reliably reconstructed barrier for NW propagation, where stress accumulates at the end of aseismic transients (Cheloni et al. 2017;Fig. 6). This barrier is the same acting on historical seismicity, hindering its general SE-to-NW migration (Napolitano et al. 2020a). The evolution of the Pollino swarm sequence, with transient aseismic and seismic events alternating west and east of the low-attenuation low-seismicity structure, suggests a mechanism of fluid migration within sealed structures developing mostly across the strike of late Pleistocene-Holocene faults (Brozzetti et al. 2017). The stress changes caused by seismic activity and the increased fault zone permeability due to fracture opening leads to fluid re-distribution through seismic and aseismic slip (Passarelli et al. 2015;Cheloni et al. 2017). The fluid distribution returns to equilibrium when the seismic activity fades out. The recovered attenuation structure represents a temporal average of the short-and long-term processes that developed in the area.
The Mercure reservoir is likely fluid-saturated so any fluid input within the reservoir increases pore fluid pressure in fluid-filled cracks around faults (Amoroso et al. 2014;Passarelli et al. 2015). The highest-magnitude seismicity consistently develops at the interface of the low-attenuation volume separating the cluster 1 and cluster 2 events, in the region of strongest short-scale cross faulting. We infer that this volume is a major structural control for the dynamics of this and following seismic sequences; if our inference about fault patches is correct, the volume could fail seismically in the future. Additional geophysical studies which integrate a range of interdisciplinary methods, a more detailed velocity model of the area, and broader seismic data sets better encompassing the study area are needed to investigate the proposed hypotheses further.

C O N C L U S I O N S
High seismic-attenuation anomalies are a reliable marker for fluid accumulation and propagation in the lithosphere. In a setting where fluid migration and changes in pore pressure are expected as the main trigger for aseismic transient and seismicity, our results offer a 3-D structural framework to interpret where fluids come from and how they interact with complex tectonics to nucleate swarms and feed springs at the surface. This fluid accumulation (depth of 7.5 km) underlies the two primary sedimentary basins of the study area, Mercure and Castrovillari.
The structures and dynamics affecting these two reservoirs appear different, with deeper high-attenuation anomalies (12.5 km) sourcing fluids only under the Castrovillari basin. There is no highattenuation volume directly below the cluster 1 events. The lowattenuation volume directly above the area of the cluster 1 events correlates with the presence of impermeable formations possibly enhanced by changes in fault mechanisms. This low-attenuation volume as well as those corresponding to the ductile low-V P eastern unit of the Apennine Platform are contoured by deep-sourced springs and likely control the spring feeding mechanisms. Seismicity before the swarm develops mostly horizontally and below the impermeable formations, supporting an interpretation in terms of lateral fluid migration. The primary deep source region for these fluids is located under the Castrovillari basin, with fluid migration possibly developing from SE (Castrovillari basin) to the NW (Mercure basin). This trend corresponds to the fluid/stress migration derived from historical seismicity and stops under the Lauria mountains due to a sharp high-to-low attenuation contrast. However, along this trend, the two clusters of the Pollino sequence develop on either side of a deep low-attenuation volume, showing no direct high-attenuation pathways for fluid migration. Our results suggest no barrier between reservoirs. More likely, delayed lateral stress exchange accommodates on small fault patches that allow fluid propagation and could break in the future.

A C K N O W L E D G E M E N T S
This work was undertaken as part of the Natural Environment Research Council (NERC) Centre for Doctoral Training (CDT) in Oil and Gas [grant number NEM00578X/1]. It is sponsored by University of Aberdeen whose support is gratefully acknowledged. We thank EIDA Data Archives for data archiving and the Universitá della Calabria for providing the data set for this study. We thank Dr Cristina Totaro for providing the velocity model used for comparison with our results. We would also like to thank two anonymous reviewers whose constructive feedback greatly improved the original manuscript.

DATA AVA I L A B I L I T Y
Waveforms of temporary stations used in this paper were provided by Universitá della Calabria and are available upon request. Seismograms recorded by INGV seismic stations (CUC, MTSN, MGR, SALB) were extracted by EIDA Data Archives at http://www.webdc.eu/webdc3/ (last access August 2017). Figures of this paper were produced using Voxler (Golden Software LLC, www.goldensoftware.com), Inkscape for MacOS X (version 1.0beta) and Adobe R Photoshop R .

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online. Figure S1. Cross-sections showing the Q −1 P variations at three frequency bands: 3, 6 and 15 Hz for Data set 1. The top single panel shows a map of the area with the seismic events used in this study (red dots) and the seismic stations (blue triangles). Line AA' represents an NW-SE cross-section (second row panels). Line BB' represents an N-S cross-section (third row panels) and line CC' represents an E-W cross-section (bottom row panels), all of them taken through the main swarm. The black circle markers represent the events of the main seismic swarm, the hourglass markers represent the events of the secondary swarm and the open square markers represent the rest of the seismicity in the area. The grey areas are beyond our resolution power and they have been masked for simplification. Figure S2. Input (right panels) and results (left panels) of two spike tests performed using Data set 1 at a frequency of 3 Hz at 7.5 km depth. The black lines on the panels represent the main faults in the area, adapted from Brozzetti et al. (2017). Figure S3. Cross-sections showing the V P (left column) and V S (right column) velocity structure in our area of interest at 3, 7.5 and 12.5 km. The model has been adapted from Totaro et al. (2014). The depth is shown at the bottom left corner of each panel. The black lines represent the main faults in the area, adapted from Brozzetti et al. (2017).
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.