Qβ, Qc, Qi, Qs of the Gargano Promontory (Southern Italy)

We have provided the first estimate of scattering and intrinsic attenuation for the Gargano Promontory (Southern Italy) analyzing 190 local earthquakes with ML ranging from 1.0 to 2.8. To separate the intrinsic Qi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{i}$$\end{document} and scattering Qs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{s}$$\end{document} quality factors with the Wennerberg approach (1993), we have measured the direct S waves and coda quality factors (Qβ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{\beta }$$\end{document}, Qc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{c}$$\end{document}) in the same volume of crust. Qβ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{\beta }$$\end{document} parameter is derived with the coda normalization method (Aki 1980) and Qc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{c}$$\end{document} factor is derived with the coda envelope decay method (Sato 1977). We selected the coda envelope by performing an automatic picking procedure from Tstart=1.5TS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{\mathrm{start}}=1.5{T}_{S}$$\end{document} up to 30 s after origin time (lapse time TL\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{L}$$\end{document}). All the obtained quality factors clearly increase with frequency. The Qc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{c}$$\end{document} values correspond to those recently obtained for the area. The estimated Qi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{i}$$\end{document} are comparable to the Qc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{c}$$\end{document} at all frequencies and range between 100 and 1000. The Qs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{s}$$\end{document} parameter shows higher values than Qi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{i}$$\end{document}, except for 8 Hz, where the two estimates are closer. This implies a predominance of intrinsic attenuation over the scattering attenuation. Furthermore, the similarity between Qi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{i}$$\end{document} and Qc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{c}$$\end{document} allows us to interpret the high Qc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{c}$$\end{document} anomaly previously found in the northern Gargano Promontory up to a depth of 24 km, as a volume of crust characterized by very low seismic dumping produced by conversion of seismic energy into heat. Moreover, most of the earthquake foci fall in high Qi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{i}$$\end{document} areas, indicating lower level of anelastic dumping and a brittle behavior of rocks.


Introduction
It is generally accepted that the seismic-wave attenuation (except geometrical spreading) consists of two mechanisms: intrinsic attenuation and scattering.The intrinsic attenuation of the direct wave is the anelastic damping due to the conversion of kinetic energy into heat.The scattering is the mechanism of elastic reflection, refraction, and diffraction of a seismic wave due to propagation through inhomogeneities or to the presence of cracks and defects in the medium.The essential difference between these two mechanisms of energy loss is that scattering attenuation is characterized by the scale length of heterogeneities (the mean free path), whereas the intrinsic attenuation is characterized by time scale (decay time) (Yoshimoto 2000).Thanks to a great intuition, Aki (1969) proposed for the first time that the seismic coda waves of local earthquakes are backscattered waves from numerous randomly distributed heterogeneities in the Earth and may be treated by statistical methods.Some years later, Aki and Chouet (1975), by using data from spectral seismograms, presented evidence for supporting Aki's idea and introduced the parameter coda Q c as a measure of the amplitude decay rate of coda waves within a given frequency band.They also showed that this decay rate was independent of recording site and event location, if observations were made within approximately 100 km from the epicenter.Starting from these pioneering papers, many studies have been carried out at inferring the attenuation of elastic waves in the lithosphere and Earth's interior by using coda waves (reference in Fehler and Sato 2003).The scattering attenuation is generally treated by three different assumptions: the single-scattering model, usually applied in the analysis of the early part of the coda (Sato 1977); the multiple-scattering model, suitable for simulating the entire coda and treated by the radiative transfer theory (Zeng 1991); and the diffusion model, which is an extreme limit of the multiple scattering process (Dainty and Toksöz 1981).In the singlescattering model, the traveling distance of the direct wave should be less than the mean free path, the scattered field is assumed weak respect to the unscattered field (Born approximation), and the source-receiver separation is neglected (Sato 1977).This approximation is then well applicable to short-period seismograms and small coda, so in this paper we will refer to the single-scattering approximation to retrieve Q c .The evaluation of the separated effects of scattering and intrinsic attenuation on coda waves is generally obtained by two approaches, both based on the multiple-scattering assumption: the multiple scattering lapse time window analysis (Fehler et al. 1992) which analyzes the coda waves over time and the Wennerberg method (Wennerberg 1993) which needs the independent measurements of the quality factors of the direct wave Q and of the coda waves Q c , sampled for the same Earth volume.In the last three decades, these two approaches were applied and developed to evaluate intrinsic attenuation and scattering in several areas worldwide both in tectonic contests (Akinci et al. 2020;Bianco et al. 2002;Londoño et al. 2022;Sharma et al. 2015;Shengelia et al. 2020;Talebi et al. 2021, among the others) and in volcanic ones (Castro-Melgar et al. 2021;Del Pezzo et al. 2019;Ibáñez et al. 2020;Prudencio et al. 2015;Ugalde et al. 2010, among the others) by estimating the quality factors Q i and Q s , respectively.
In this paper, the separation of Q s and Q i is attempted for the Gargano Promontory (Southern Italy, hereafter GP).The Wennerberg approach has been used, and this implies the independent estimations of Q and of Q c for the same area.Q has been estimated using the coda normalization method (Aki 1980) and Q c has been estimated using the method of the coda envelope decay (Sato 1977), over the same dataset as Q .The GP is classified as high seismic hazard (horizontal peak ground acceleration PGA = [0.200− 0.225] g for the 10% probability of exceedance in 50 years, from MPS04 of Stucchi et al. 2004).This classification follows the 2002 October 31th earthquake ( M W = 5.7 ) that hit Molise (South- ern Italy), a region bordering the GP and until then classified as non-seismic.In 2013 the GP was covered by the OTRIONS local seismic network, installed and managed in cooperation between the University of Bari Aldo Moro (hereafter UniBa) and the National Institute of Volcanology and Geophysics (hereafter INGV), and since then an intense microseismicity has been recorded.The present work benefited of the recently released database of waveforms of the GP microearthquakes (Filippucci et al. 2021c).

Geological and seismological setting of the studied area
From a geodynamic point of view, the GP is considered a continental lithospheric block belonging to Adria plate, subducting the Apennines chain in the western side (de Lorenzo et al. 2017;Petrullo et al. 2017;Bentivenga et al. 2017;Bucci et al. 2019) and Dinarids and Albanids in the eastern sector (Di Bucci and Angeloni 2013).
The structural setting of GP (Fig. 1) is complex and subject to different interpretations.Ortolani and Pagliuca (1987)  The GP area is adjacent to the Apennine chain which is the most seismically active part of the Italian territory, both for frequency recurrence and magnitude of earthquakes.Even if the seismic history of GP is characterized by destructive events, as the 1627 earthquake (estimated M W = 6.7 ) and the 1646 earthquake (estimated M W = 6.2 ) (Del Gaudio et al. 2007), in the last 5 decades, few events with magnitude M W > 4 have been recorded in the GP.The number of high-magnitude earthquakes increases if we move some tenths of kilometers westward toward the Molise, where the 2002 October 31th earthquake occurred (Rovida et al., 2020;Rovida et al. 2022).
Since the installation of the first seismic station on the GP (by INGV and UniBa in 1986), the instrumental coverage of the GP has been improved up to 4 stations in 2008.In 2013, UniBa and INGV cooperated for the installation of 12 new seismic stations (OT network) which greatly improved the seismic Fig. 1 Seismotectonic and geological maps of the studied area.Left: Adria microplate (simplified from Mantovani et al. (2006).The GP is enclosed in the dashed line square.Right: geological map of GP, modified by Miccolis et al. (2021).The colors of the stratigraphic units are related to the geological ages (Jurassic, light blue; Cretaceous, green; Eocene, orange; Quaternary, white).The igneous outcrop of Punta Pietre Nere, cited in the text, is marked with a black star.The acronyms refer to the main faults (MF, Mattinata Fault; AF, Apricena Fault; CF, Candelaro Fault; SF, Sannicandro Fault) Vol.: (0123456789) monitoring of the southern Italy and of GP.In the following years the network has undergone some changes (refer to Filippucci et al. (2021c) for more details) up to the current layout, which can be viewed on the website eida.ingv.it/it/networks/network/OT(University of Bari "Aldo Moro" 2013).The seismic monitoring of OT and IV networks on GP revealed an intense activity of micro-earthquakes with magnitude up to 2.8 and the results are disclosed in literature by some recent papers.The 1D velocity model (de Lorenzo et al. 2017) revealed that the GP seismic activity is concentrated at two different depths: deeper in the NE with foci below 20 km and shallower in the SW with foci above 15 km.This difference is probably due to a slope of the brittle/ductile transition at depth, as initially hypothesized by Filippucci et al. (2019b) and subsequently numerically modeled by Lavecchia et al. (2022).The hypocenter distribution, moving from southwest to northeast, aligns itself by drawing a seismogenic surface that deepens toward the northeast but on which earthquakes rupture in an orthogonal direction, toward the northwest, with thrust faulting mechanism (Filippucci et al., 2020;Miccolis et al. 2021).The depth of this lineament agrees with the Moho depth in the GP which is higher (25-30 km as estimated by de Lorenzo et al. (2014;2017)) compared to the Apennine (18-20 km) (Cassinis et al. 2007).The active faults responsible of the seismicity of GP are not yet recognized and no relationship has been found with the San Marco in Lamis-Mattinata composite fault reported by the DISS 3.3.0database (DISS Working Group 2021).The southernmost part of GP is characterized by fluid circulation along fractures at shallow depth (Tripaldi 2020), and this may be the origin of the shallower seismicity of this area.In this frame, Filippucci et al. (2019a), by carrying out a 2D attenuation study at a local scale, revealed that along the direction in which the seismicity becomes deeper, that is moving from southwest to northeast, Q c decreases.These authors found also a high Q c anomaly which can be correlated with the anomalies of other geophysical observations as gravity (Loddo et al. 1996), surface heat flow (Vedova et al. 2001), V P ∕V s (Improta et al. 2014), and mag- netic field (Speranza and Chiappini 2002), but they were not able to interpret it in terms of anelasticity or scattering.In order to answer to the question whether these anomalies are related to anelasticity or to heterogeneities of the GP crust, in this paper we separate the contribution of quality factors Q i and Q s to coda attenuation and finally discuss the trend of Q c with depth recently obtained by Filippucci et al. (2021b).

Earthquake dataset
The dataset is shown in Fig. 2 and it consists of 190 low magnitude earthquakes (1 ≤ M L ≤ 2.8, histogram in Fig. 2D) recorded between July 2015 and August 2018 by the UniBa OTRIONS (FDSN code OT) and the INGV RSN (FDSN code IV) seismic networks (refer to Filippucci et al. (2021c) for a complete description of the OT seismic network) (Fig. 2) for a total of 24 stations.The OT network started to operate in 2013, and from 2019, recordings are transmitted continuously and in real time on the EIDA platform (European Integrated Data Archive).The OT stations are equipped with 3-component short-period Lennartz velocimeters (Filippucci et al. 2021c); the IV stations are equipped with 3-component velocimeters ranging from very broadband to short period (Michelini et al. 2016).All the considered earthquakes are localized in the GP, with epicentral distance less than 70 km, foci depth reaching the lower crust down to 35 km, and M L ranging from 1.0 to 2.8 (Fig. 2).Locations have an average horizontal error of 0.85 km and an average vertical error of 0.94 km.Events with hypocentral distances larger than 70 km were removed from the analysis being out of the area of interest.This dataset has been already released and is available for downloading (Filippucci et al. 2021d).
The data set of earthquakes used in this paper was already presented in a previous paper (Filippucci et al.Gaudio et al. 2007) included in the dotted boxes (see text for more details).B Record section (Palombo et al. 2022) with the ~ 30% of the analyzed waveforms (one every 0.2 km).Each one is normalized by its absolute maximum amplitude.Seismograms are sorted by hypocentral distances (y-axis) and aligned by origin time (at 0 s in the x-axis).The green curve denotes direct S-wave arrival times (T S ), the light blue curve indicates 1.5 T S , and the red one shows the lapse time (T L ) at 30 s from origin time.C Percentage histogram of earthquake magnitudes.D Percentage histograms of earthquake foci depths ◂ Vol:.(1234567890) 2021b).In this article, we selected a dataset of waveforms and envelopes that was the result of a manual selection, which was very time expensive.The SNR was a posteriori computed and we verified that for all the waveform SNR > 2 in all frequency bands.Moreover, all the waveforms and the envelopes have been published and released (for waveforms : Filippucci et al. 2021c;Filippucci et al. 2021d;for envelopes: Filippucci et al. 2021e;Filippucci et al. 2021a).
A total of 1477 3-component recordings were initially considered.Then, we developed an automatic procedure for waveform picking to accept all the recordings for which T L > 1.5T S , being T L the time measured from the earthquake origin time and T S the S wave travel time.
In Fig. 2, the record section on the right shows that the choice of T L = 30s allows to include in the dataset the seismograms recorded at hypocentral distance less than 70 km.For T L = 30s , the radius of the spherical volume of Earth crust involved in the Q c computa- tion is 58 km for GP (as already derived by Filippucci et al. (2019a)).With T L = 30s , a total number of 1317 3-component recordings were suitable for further analysis, that is, the 89% of the initial dataset.

Q β estimation
We estimated Q applying the coda normalization method (Aki 1980).This method adopts the ratio between direct S wave and coda wave amplitude spectra at a fixed lapse time T L .The obtained spectral ratio is considered independent of site and S-wave source effects (Sato et al. 2012).
Q was estimated for several frequencies f = {3, 4, 5, 6, 8, 10, 12} Hz with frequency bands ranging from f ∕ √ 2 to √ 2f by performing a linear regression of log- arithmic normalized amplitude vs. distance, based on the following relation (Aki 1980): where r is the hypocentral distance, A β and A c are, respectively, the spectral amplitudes of S wave and coda waves, Δr = 1.25 km, V = 3.86 km/s is the S wave average velocity in the GP (de Lorenzo et al. (1) ⟨ln 2017), and a(f ) is the intercept of the linear fitting.Geometrical spreading is usually represented in the left side of Eq. (1) as r , being the geometrical spreading factor which varies depending on the considered seismic wave and source-receiver distance (Frankel 2015 and bibliography therein).For body waves and focal depths up to 100 km, = 1 (Yoshi- moto et al. 1993; Frankel 2015).
We set T L = 30 s and fixed the time window length Δt w = 2.56 s.Δt w starts 0.1 s before T S for measuring A ; Δt w is centered on the lapse time T L , over coda waves for measuring A c .
This choice of the time window length causes the number of sampling points to be a power of 2, being the sampling rate of the OT network seismometers at 100 Hz.The discrete Fourier transform in the spectral analysis module of internal functions of SAC (Goldstein and Snoke 2005; Goldstein et al. 2003) was computed without zero padding to obtain the S-wave and coda wave spectra.A 5% cosine taper was applied to signals inside Δt w in order to reduce dis- tortion caused by finite length signals and avoid the Gibb's phenomenon (Scherbaum 2001).Q estimates were performed without using different time window amplitudes since Yoshimoto et al. (1993) showed that Q values vary slightly with the length of the time win- dows, always within the standard deviation.
We considered the three components of the seismograms filtered with 4-pole Butterworth filter in the frequency band defined by � .So, the considered frequency bands (in Hz) are f = 3 Hz: [2.1-4.2],f = 4 Hz: [2.8-5.7],f = 5 Hz [3.5-7.1],f = 6 Hz [4.2-8.5],f = 8 Hz [5.7-11.3],f = 10 Hz [7.1-14.1],f = 12 Hz [8.5-17].In Fig. 3, the values of ln rA (r,f ) , grouped in space windows of Δr = 2.5 km, are plotted as function of r for each frequency band.Several tests were carried out to obtain the optimal choice of ∆r.Since Δr has an effect of smoothing on data, we inferred that Δr = 2.5 km is the best com- promise between the need of obtaining an acceptable value of the determination coefficient ( R 2 > 0.7) and the need of minimizing the error on Q .Q (f ) is the result of the linear regression through Eq. ( 1), and it is plotted as a function of frequency in Fig. 4 and listed in Table 1 with the relative errors ΔQ .
To analyze the robustness of the linear fitting in Fig. 3, we compared the results of Q obtained  1) vs. the source to receiver hypocentral distance r at each frequency.The Q β parameters, obtained by the linear fitting, are reported with their error.All the fitting parameters are reported in Table 1 Vol:.( 1234567890) by using Eq. ( 1) with those obtained by quadratic approximation of the coda normalization method (de Lorenzo et al. 2013b) arising from the inversion of log normalized amplitudes.As described in de Lorenzo et al. (2013b), the preferred Q value, at each frequency, can be obtained by comparing the values of the Akaike AIC statistical parameter: where N t is the number of data points, for which the misfit line is less than the error on data, N s is the num- ber of data points (always 28 as showed in Fig. 3), N m is the number of unknown model parameters ( N m = 2 for linear fitting, N m = 3 for quadratic fitting), and E is the variance.
Repeating the test for all the frequencies (Table 2), we obtain always that the AIC for the linear fitting is less than the AIC for the quadratic one, demonstrating that the linear fit is always more reliable than the quadratic fit.The Akaike test indicates that the increase in model complexity, in this case the quadratic approximation, does not imply a better fit of data.

Processing for Q c estimation
We performed a new estimation of Q c by consider- ing the same source to receiver couples previously used for inferring Q .We applied the method of coda envelope decay (Aki 1969) based on the linear fitting of coda energy envelope.An automatic procedure was implemented to pick the starting time of envelope decay of coda waves.The procedure consists of selecting coda waves in a time window between 1.5T S and T L , being T L = 30 s.In our case, the choice of 1.5T S allows to include recordings at distances less than 70 km (Fig. 2).It is well known that the estimated Q c is influenced by the time win- dow used for linear regression.In fact, depending on this choice, coda waves sample different regions of the lithosphere.The early coda corresponds to a  smaller sampling ellipsoid (centered at the hypocenter and the station).This is the reason for the discrepancy between the results of our study and the results of the preceding study (Filippucci et al. 2021a, b, c, d, e) as shown in Fig. 5.We used the same dataset selected for the Q esti- mate and observed that all the 1317 three-component seismograms were suitable for the Q c analysis.The waveforms were bandpass filtered with fourth-order Butterworth around the central frequency f = {3, 4, 5, 6, 8, 10, 12} Hz, with frequency bands defined in the previous subsection.Assuming the single isotropic scattering theory (Sato 1977), Q c (f ) can be retrieved by using the following relationship (Sato 1977): (2) ln 2) allows to estimate Q c (f ) from a linear regression analysis of the decay rate of the logarithmic coda envelope A c (r, f , t) vs. time t .Outliers generally correspond to small values of the coefficient of determination ( R 2 < 0.01 ) and have been removed from the dataset.
Usually (e.g., Wang and Shearer 2019) the coefficient of correlation R is computed on smoothed enve- lopes, and for this reason a R > 0.9 R 2 > 0.81 is achievable.In our study, we used R 2 > 0.01 ( R > 0.1 ) as criterion of data removal, but on raw envelopes.In fact, at any frequency, the envelope smoothing procedure does not affect Q c , whereas it can drastically increase R 2 , until it reaches a statistically acceptable  Vol:.( 1234567890) value.Anyway, the smoothing procedure, eliminating the casual fluctuation of the data, can produce a lower variance on Q c .Therefore, we preferred to work on raw envelopes.
The averaged estimates of Q c (f ) are plotted in Fig. 4 and listed in Table 1 with the standard deviation ΔQ c .

Separation of Q i and Q s effects
We applied the method proposed by Wennerberg (1993) to infer the scattering and the intrinsic quality factors, Q s and Q i , separately.The Wennerberg's (1993) approach is valid for the case of co-located source and receiver; it is based on the intuition that the ratio between the coda amplitude decay predicted by the multiple scattering model (Zeng 1991) and that predicted by the single scattering model (Aki and Chouet 1975) follows a linear trend as a function of the dimensionless mean free time = T L ∕Q s , where = 2 f .Assuming that Q c is computed for the single- scattering model, Q s and Q i can be derived for the multiple-scattering model by this mathematical relationship: being ( ) a correction needed to account for multi- ple-scattering model of Zeng (1991).( ) is a theoreti- cal discrete point function of that can be approximated by a linear function (Wennerberg 1993), where ( ) is given by (Del Pezzo et al. 1995)   According with Wennerberg (1993), if we assume that Q is the total attenuation Q T , the combina- tion of intrinsic and scattering attenuation effects is described by The necessary condition to separate the contributions of intrinsic and scattering attenuation is to have independent estimates of Q c and Q , both obtained as function of frequency for the same Earth volume. (3) This is the reason why we computed Q c instead of using previous results (Filippucci et al. 2019a(Filippucci et al. , 2021a)).Then, one approach is to solve the system formed by Eqs. ( 3) and ( 5), which leads to Since Q T = Q we can also evaluate the seismic albedo T defined as the dimensionless ratio of scattering loss to total attenuation, and the extinction length L e = v 2 f Q T that is the distance over which the primary S wave energy is decreased by e −1 .
The obtained results and relative errors for T L = 30 s are reported in Table 1 and plotted in Fig. 4 as functions of f .
As already discussed, the Wennerberg approach represents an approximate solution of the radiative transfer equation based on the single-scattering assumption.We checked the applicability of the Wennerberg's approximation to the present data set by quantifying the deviation of the approximate Wennerberg solution of the energy equation from the general solution of Paasschens (1997) by calculating the average percent variation (refer to apv in the Appendix A of de Lorenzo et al. ( 2013)).To compute apv we used the values of Q s and Q i , inferred in this study, at all frequencies, at distances of 30 km and 60 km, for T L = 30 s. Results indicate that apv ranges between 10 and 20% , inside the percentage error bar on Q s and Q i (Table 1), as already demonstrated for the Umbria- Marche (Central Italy) by de Lorenzo et al. (2013b).Further studies based on MLTWA may give more robust estimates of attenuation than those based on the current single-scattering hypothesis.

Dependence of Q on f
The frequency dependence of Q c , Q , Q i , Q s at T L = 30 can be retrieved by using the average value of each quality factor plotted in Fig. 4 by means of the power law relationship of rock quality parameter ( 6) We obtained the following values of Q 0 and : Q s (f ) increases with frequency with = 1.52 faster than the other quality factors, indicating that the more pronounced increase of Q s may be caused by the pre- sent half-space assumption, as already observed by Akinci et al. (2020) for the central Apennine.

Discussion
We obtained the first estimation of Q , Q i , Q s quality factors in the GP and a new estimation of Q c .We used the micro-seismicity of GP recorded by the OT and the RSN seismic network in the period 2013-2018, recently released in a database of seismic waveforms (Filippucci et al. 2021c,d).The new estimation of Q c was achieved by the implementation of an automatic procedure to pick the starting time of envelope decay over coda waves with Eq. ( 2).Since the slope of the envelope decay can vary over time (in general, the slope decreases with time, Filippucci et al. 2021a), the choice of the starting point T start could control the Q c results.Theoretically, Eq. ( 2) could be applied to a time window that starts immediately after T S but this is not the common choice.Commonly, T start = 2T S (e.g.Jin and Aki 1988;Ibáñez et al. 1990;Eyidoğan et al. 1996;Mamada et al. 1997) or in other cases T start = 1.5T S (e.g., Padhy et al. 2011 and references therein).Mukhopadhyay and Sharma (2010), for the Himalayan region, analyzed the effect of T start on attenuation and observed non-significative differences on Q c for different T start ; significative decrease in the total number of available waveforms as T start increases from 1.5T S to 2.5T S , especially for lapse times below 50 s.So, T start = 1.5T S is the most conservative choice in terms of data selection.To avoid "a priori" choices on T start and to eliminate envelopes with abrupt vari- ations along the coda decay (e.g., signal bumps), some authors prefer to manually pick the decay onset by visually inspecting the envelopes (de Lorenzo et al. 2013b;Filippucci et al. 2019aFilippucci et al. , 2021a)).The ( 8) comparison between Q c retrieved from the same data- set but with a manual picking of T start (Filippucci et al.  2021a) and Q c of this work reveals that T start = 1.5T S systematically falls into a zone of the envelope with a less decay rate respect to the manual picking and produces a small shift between the two estimates inside the standard deviation.This effect can be observed in Fig. 5, revealing a slightly less attenuating crust for GP than the previously estimated one.
Q c is confirmed to be frequency dependent and increases with increasing frequency, as already observed worldwide both in tectonic and in volcanic contests.
Q was estimated by the coda-normalization method (Aki 1980).This method eliminates source, site, and geometrical spreading effects by normalizing S-wave amplitudes with coda waves.It reveals the attenuation on body waves due to the medium anelasticity and to elastic scattering on heterogeneities.Generally, the normalized amplitudes of S-waves decrease with increasing of hypocentral distances, demonstrating that body waves attenuate in amplitude with increasing travel distance and this attenuation depends on both frequency and geometrical spreading.This effect is shown in Fig. 3 where, from the slope of the linear fit line [Eq.( 1)], Q can be computed for different frequency bands.As it can be observed from Table 1 and Fig. 4, Q regu- larly increases in all the frequency range except for 6Hz < f < 8 Hz where it is less steep.This observa- tion can be inferred also from plots in Fig. 3.The revealed S-wave attenuation at f c = 8 Hz might be due to the scattering loss at random heterogeneities having a characteristic length or alternatively the 8 Hz frequency might be the predominant frequency band of the anelastic response of the medium, as discussed by Sil et al. (2022) for the Kathmandu region, Nepal.The S-wave attenuation was computed by the assumption of the half-space approximation with uniform V , that is quite unrealistic especially for the shallower layers in the upper crust.This simplification may produce a bias in the estimation of Q .Some authors (Akinci et al. 2020;Pisconti et al. 2015) estimated the effect of a depth-dependent crustal model on the attenuation for the Central Apennine (Italy).They found out that Q s and Q i are slightly affected by these assumptions (depth-dependent or homogeneous Earth) in the case of a continental crust, where Vol:.( 1234567890) shallow earthquakes occur, while, with increasing the source-receiver distance, the observed apparent attenuation decreases, suggesting an increasing of propagation efficiency with depth below the Moho.The seismicity of GP occurs at crustal depth above the Moho discontinuity (Miccolis et al. 2021) Q attenua- tion should not be severely affected by the homogeneous half-space approximation.

Comparison with Q s and Q i estimates in other areas
The separated contribution on coda waves of anelastic attenuation Q i and elastic scattering Q s was inferred by using the Wennerberg (1993) method.The application of this method requires the availability of independent estimates of Q and Q c as functions of fre- quency for the same Earth volume in order to avoid bias in the results.In our case we used the same dataset, representative of the same crust volume of the GP, for the estimates of the total attenuation of body waves Q and of the coda wave attenuation Q c in the single-scattering approximation.Results in Fig. 4 indicate a regular ascending trend of Q c and Q i as functions of frequency.Q s presents an irregular low value at 8 Hz reinforcing the idea that the lower value of Q at 8 Hz might be due to a scattering effect.Gen- erally, Q c regularly increases with frequency.In the assumptions that the autocorrelation function of the elastic heterogeneities in the Earth medium follows a Gaussian (or von Karman) statistics, their power spectral density function shows a regular pattern with frequency, implicitly justifying a regular pattern of Q c as a function of frequency (for a wide and thorough discussion see Sato et al. (2012)).An anomaly in the heterogeneous size distribution (as compared with a "regular" distribution, Gaussian, or von Karman) may produce the observed frequency anomaly in Q s pattern.
Assuming a uniform half space with constant velocity V = 3.86 km/s as we did, f = 8 Hz corre- sponds to a wavelength ≈ 0.5 km, which could repre- sent the average distance of the random scatterers in the GP.
For all the other frequencies, scattering on heterogeneities turns out to be less important than anelasticity in the coda attenuation in agreement with the preceding results in regional studies of Southern Apennine which included marginally the GP area (Bianco et al. 2002).
Results of Q i and Q s in tectonic contests worldwide suggest that the attenuation of coda waves is dominated by the intrinsic attenuation (Padhy et al. 2011;Farrokhi et al. 2016;Pujades et al. 1997;Akinci et al. 1995;Aki and Chouet 1975;Bianco et al. 2005;Tuvè et al. 2006, among others).The predominance of the anelasticity on coda waves increases by increasing frequency since the percentage ratio (Q s − Q i )∕Q s moves from 48 at 3 Hz to 70% at 12 Hz (Table 1).The seismic albedo is B 0 < 0.5 for all the frequency bands, which indicates that anelasticity is the predominant attenuation effect in the region.The only exception is for f = 8 Hz where B 0 ≅ 0.5 , indicat- ing a predominance of scattering at the scale length of this frequency.The extinction length L e ranges between 18 and 31 km, increases with frequency (with the exception for f = 8 Hz where L e is lower), and it is comparable with other studies in tectonic domains (Akinci et al. 2020;Londoño et al. 2022;Sharma et al. 2015;Shengelia et al. 2020;Talebi et al. 2021, among others).
The higher is T L , the longer is coda and the num- ber of waveforms, as can be observed in Fig. 2B.In our study, the choice of T L = 30 s is a compromise between the need of high number of coda signals with short duration.This choice implies that the selected coda signal can include contributions of S waves traveling in the mantle.This occurs because if we consider an average depth of hypocenters about 20 km and a time window of T L = 30 s, S coda waves, trave- ling with V = 3.86 km/s, can arrive from a maximum depth of about 68 km, in the case of vertical incidence for sake of simplicity.Considering a Moho depth of about 40 km, our results might be biased by the constant velocity half space assumption (Del Pezzo et al. 2019).A wide discussions on this topic can be found in Margerin et al. (1998Margerin et al. ( , 1999)), Del Pezzo and Bianco (2010), and references therein.These authors demonstrated that if the scattering strengths of the mantle and the crust are comparable, the presence of a velocity contrast at the Moho could amplify the coda signal, since part of the energy remains trapped in the crust; if the scattering strength of the mantle is smaller than that of the crust, the shape of the envelope decay can be reproduced by the Wennerberg (1993) model (Margerin et al. 1998).It is also worth to note that a scattering strength for the upper mantle smaller than that of the crust was previously inferred for the Italian peninsula (Bianco et al. 2005).Further studies based on a more realistic velocity model (e.g., del Pezzo et al. 2011) can be used to discriminate the scattering strength of the crust from that of the mantle.In the case of the GP area, the Moho depth is not well constrained.The Moho is thought to be located between 35 and 40 km (Chiarabba et al. 2005;Piana Agostinetti and Amato 2009) or at depths greater than 40 km (Mele and Sandvol 2003).In fact, at depth greater than 30 km, the seismicity become sparse (less than 5%, Fig. 2D), the velocity contrast due to the Moho is not observed (de Lorenzo et al. 2017), and the trend of seismicity, that deepens toward the Adriatic Sea, probably follows the trend of the deepening of the Adriatic Moho (Miccolis et al. 2021).
The comparison between the outcomes of this study and those previously obtained for other tectonically active regions in Italy and worldwide reveals that the GP values of seismic quality factors Q i and Q s are within the end members at all frequencies (Fig. 6).

New insights concerning the geophysical anomalies of the GP
The frequency dependence relationships of Q , Q c , Q i , Q s show estimates of Q 0 up to 40 and val- ues of greater than 1.25.The paradox consists of both the small values of Q 0 (less than 48) and the high values of (greater than 1.1) that should indicate a high tectonic activity and the general absence of high-magnitude earthquakes in the GP (Morozov 2008).The M W = 5.7 Molise earthquake (Southern Italy, 2002 October 31th) occurred in an adjacent area previously classified as low hazard area, may reinforce the hypothesis that the GP could reasonably located in a tectonically active regime.This hypothesis is also sustained by the observation that more than ten historical earthquakes with estimated magnitude Mw > 5.5 occurred in the GP and its surrounding in the last millennium (Fig. 1A).
At 12 Hz we can observe in Fig. 4 the maximum difference between Q s and Q i , so it might be stated that coda amplitudes at 12 Hz are controlled by anelastic damping and it can be considered predominantly attenuation of body waves, since attenuation of surface waves is detectable at lower frequencies, less than 10 Hz (Aki and Chouet 1975).Therefore, Q c ≈ Q i at all frequencies and this result has relevant consequences if we examine the recent 3D tomography of coda attenuation in GP (Filippucci et al. 2021a).The 3D images of Q c were obtained for three frequencies (3 Hz, 6 Hz, and 12 Hz in Filippucci et al. (2021a)) by using the Del Pezzo and Ibáñez (2020) approach which consists of computing the polynomial approximation of the analytical sensitivity kernels.The polynomial sensitivity kernels were then used as weighting functions for estimating Q c in an Earth vol- ume, which was subdivided in cubic pixels with side of 5 km, producing a horizontal Q c image every 8 km of depth for the GP crust.We reported in Fig. 7 the images of the 3D tomography of Q c at 12 Hz (modi- fied from Filippucci et al. 2021a).Now we can consider Q c ≈ Q i and give a physical interpretation of the observed anomalies by interpreting them as attenuation of body waves.The well-defined high Q i anom- aly, extending down to 24 km, reveals the presence of a body embedded in a more attenuating one in the northern GP which might agree with the hypothesis of the existence of a high-density and high-susceptibility body in the same area (Loddo et al. 1996).This hypothesis agrees also with the observation of some igneous rocks that crop out in GP in a land named Punta delle Pietre Nere (black star in Fig. 1), intruding the upper Triassic sedimentary successions.These magmas are interpreted as derived from an amphibole-bearing lithospheric mantle source at 70-90 km depth (Mazzeo et al. 2018).
Figure 7 shows a Q i decreasing trend (attenuation increasing trend) moving toward the northeast sector of GP.This sector is characterized by an anomalous absence of seismic activity down to 20 km (Miccolis et al. 2021).Therefore, the observed Q i decreasing trend well correlates with the hypothesis of a ductile behavior characterizing the upper basement and the sedimentary cover, as confirmed by a recent thermorheological model (Lavecchia et al. 2022).At depths greater than 20 km, a seismogenic layer is encountered (Miccolis et al. 2021) and it was modeled as due to the presence of fluids in the deepest part of the basement (Lavecchia et al. 2022).These results agree with the presence of general high Q i (low attenuation) values that can be a hint of a brittle and low-strength lower crust in the deepest seismogenic layer (Fig. 7, 32 km).The correlation between seismic attenuation and rock rheology has been widely explored by several authors by using numerical methods (Castillo et al. 2022;Castaldo et al. 2019, among others) and laboratory experiments (Sato et al. 1988;Carcione et al. 2020, and references therein).
It is worth to note that at all depths, the most part of earthquake foci falls in areas characterized by high Q i , indicating a quite good agreement between the presence of seismicity and lower levels of anelastic damping both correlated to the brittle behavior of rocks.

Conclusions
In the list below, we summarize the results of this study and the main implications and observations already detailed in the discussion section: • The installation of the OTRIONS seismic network and its inclusion in the national array managed by INGV allowed us to record many low-magnitude earthquakes in the GP area in the last decade.• The first estimate of the attenuation of shear waves has therefore been obtained for the GP using the coda normalization method.The linear fit is robust, as shown by comparing the obtained results with those arising from a quadratic approximation of the logarithmic ratios of amplitudes.• After the study of Filippucci et al. (2021b), Q c has been newly estimated by assuming as starting time of the coda window T start = 1.5T S and results indicate that the value of this parameter depends on the used T start , as discussed in several previous studies.• The first estimate of scattering Q s and intrinsic Q i quality factors has been obtained for the GP.The main assumptions (single scattering, homogeneous half-space) of our approach may cause a bias on the obtained Q s and Q i values.Comparison with the solution obtained with the multiple-scattering hypothesis reveals that these biases should be of the same order as the uncertainties associated with the current estimates of Q .However, thanks to the improvements of the OTRIONS seismic network and the ever increasing number of good quality data, further studies including the multiple scattering hypothesis can be conducted with the dual purpose of evaluating its effect on Q s and Q i and of calculat- ing the scattering strength of the crust and mantle.• The low value of Q s at f = 8 Hz can be an indica- tion of the presence of random scatterers in the GP crust with an average distance of ≈ 0.5 km.• The comparison of the Q estimates in the GP area with those obtained in the Umbria-Marche (de Lorenzo et al. 2013b), under the same assumptions, reveals that the elastic and anelastic properties of the GP should be somewhat different from those of adjacent Apennines, indicating a different stressstrain regime that may result in a different level of seismicity of the two tectonically active areas.• Our results indicate that Q c ∼ Q i , so we are now able to obtain a physical interpretation in terms of inelastic attenuation of the magnetic, gravimetric, and coda attenuation anomalies found in previous studies in GP.
consider central portion of GP a pushup structure influenced by Dinarids.Patacca et al. (2008) invoke halokinetic factors along with tectonics.Bertotti et al. (1999) propose a structural model up to a depth of 4 km, based on surface geological evidence, seismic, and well data.These authors hypothesize a basal detachment of the entire carbonate multilayer, located at a depth of 4 km in the Carpino location and becoming more superficial moving toward SW.Near Monte Granata the detachment intersects the topographic surface.

Fig. 2 A
Fig. 2 A Earthquakes recorded by OTRIONS and INGV seismic networks between 2015 and 2018, processed in this study.Left: geographical map of the studied area where gray dots are earthquakes, yellow triangles are INGV stations and orange triangles are OTRIONS stations.Blue lines are source-station paths.The gray box shows the most important historic and recent earthquakes occurred in the area (from DelGaudio et al. 2007) included in the dotted boxes (see text for more details).B Record section(Palombo et al. 2022) with the ~ 30% of the analyzed waveforms (one every 0.2 km).Each one is normalized by its absolute maximum amplitude.Seismograms are sorted by hypocentral distances (y-axis) and aligned by origin time (at 0 s in the x-axis).The green curve denotes direct S-wave arrival times (T S ), the light blue curve indicates 1.5 T S , and the red one shows the lapse time (T L ) at 30 s from origin time.C Percentage histogram of earthquake magnitudes.D Percentage histograms of earthquake foci depths

Fig. 3
Fig.3Average values of the left member in Eq. (1) vs. the source to receiver hypocentral distance r at each frequency.The Q β parameters, obtained by the linear fitting, are reported with their error.All the fitting parameters are reported in Table1

Fig. 4
Fig. 4 Plot of Q , Q c , Q i , Q s quality factors versus frequency f obtained for the GP

Fig. 5
Fig. 5 Comparison between Q c parameter obtained in this study with that estimated by Filippucci et al. (2021a)

Fig. 6
Fig. 6 Comparison of a Q i and (b) Q s estimated in this study with those obtained for other tectonically active regions of Italy and worldwide.Color and style of curves are related to the method adopted by the authors to separate Q i and Q s [Wennerberg's (1993) method with gray solid line; MLTWA method (Fehler et al. 1992) with black and dashed line].The plotted