3D-Kernel Based Imaging of an Improved Estimation of (Q_c) in the Northern Apulia (Southern Italy)

Abstract

We investigate crustal seismic attenuation by the coda quality parameter ($Q_c$) in the Gargano area (Southern Italy), using a recently released dataset composed of 191 small earthquakes (1.0 ≤ ML ≤ 2.8) recorded by the local OTRIONS and the Italian INGV seismic networks, over three years of seismic monitoring. Following the single back-scattering theoretical assumption, $Q_c$ was computed using different frequencies (in the range of 2–16 Hz) and different lapse times (from 10 to 40 s). The trend of $Q_c$ vs. frequency is the same as that observed in the adjacent Umbria-Marche region. $Q_c$ at 1 Hz varies between 11 and 63, indicating that the area is characterized by active tectonics, despite the absence of high-magnitude earthquakes in recent decades. The 3D mapping procedure, based on sensitivity kernels, revealed that the Gargano Promontory is characterized by very low and homogeneous $Q_c$ at low frequencies, and by high and heterogeneous $Q_c$ at high frequencies. The lateral variations of $Q_c$ at 12 Hz follow the trend of the Moho in this region and are in good agreement with other geophysical observations.

1. Introduction

Coda waves recorded from small local earthquakes are interpreted as body waves of backscattering from heterogeneities that are distributed randomly but uniformly in the Earth’s crust and upper mantle. Coda waves are predominantly composed of scattered S-waves and, in order to relate coda wave characteristics to seismic attenuation, a scattering model is necessary. In the early 1970s, Aki and Chouet [1] and Sato [2] developed a single-scattering model in the Born approximation, where the primary wave encounters a scatterer (inhomogeneity) only once, producing a secondary wave reaching the receiver with no other impact. This model parametrizes the seismic attenuation using a single parameter, $Q_c$, called Q-coda, which is easily calculated once the decay rate of the coda envelope is measured from data. Worldwide, $Q_c$, estimates generally report low values in regions with active tectonics, and high values in stable regions [3,4,5,6].

The physical interpretation of $Q_c$ is a long-standing problem in seismology. Years after the development of the early single-scattering models, it was demonstrated by Zeng et al. [7] that the single-scattering model is a first-order approximation of the general solution to the energy transport equation (ET) ([8], which accounts for multiple scattering and is parametrized in terms of two parameters, intrinsic $Q_i$ and scattering $Q_s$. New interpretations of coda wave envelopes, based on the Paasschens solution of the ET equation, then demonstrated that single-scattering models are unrealistic for most of the experiments carried out through the Earth (see, e.g., [9,10,11] for a review), while the other asymptotic approximations of the ET equation, the so-called diffusion model [12], are much more realistic in cases showing high geological heterogeneity, like volcanoes. Despite this evident limitation, $Q_c$, or, which is the same in practice, the coda envelope decay parameter, is still used, due to the easy application of the single-scattering model to experimental data.

In the single-scattering model, $Q_c$ is a combination of scattering $Q_s$ and intrinsic $Q_i$ attenuation ($Q_c^{-1}=Q_s^{-1}+Q_i^{-1}$) [1]. It is interesting to note that the diffusion model predicts that $Q_c$ asymptotically approaches the intrinsic attenuation $Q_i$ at increasing times in the coda wave train (assuming a uniform half-space [12]). However, single-source-receiver couple estimates of $Q_i$ and $Q_s$ can only be obtained in strongly heterogeneous media, where the diffusion model is applicable. For all other cases, the separation is almost impossible and the only measurable (unbiased) parameter is $Q_c$ [11].

Sato’s [2] derivation of the single-isotropic-scattering model, and the numerical computations of Frankel and Wennerberg [13], suggest that $Q_c$ is an intrinsic attenuation parameter. Zeng’s [7] simulations imply that $Q_s^{-1}$ must be approximately from 5 to 15 times as $Q_i^{-1}$ to have a comparable effect on the coda decay rate.

To achieve a correct interpretation of the measured $Q_c$, the association of the attenuation within volumes crossed by the coda radiation with the structural, tectonic, or geodynamic maps must be performed. Starting from the early and approximative $Q_c$ images [14], many attempts were made to spatially confine the attenuation anomalies measured through coda waves (as reviewed by del Pezzo et al. [11]), in order to reveal heterogeneous Earth structures and interpret them in terms of rock quality. Among the imaging methods, sensitivity kernels are essential to obtain an attenuation image from coda waves; the theoretical approach, in the framework of both radiative transfer theory and the diffusion model, focuses on the spatial changes of $Q_c$ from distributed sources recorded at a seismic network [15,16]. In recent years, the sensitivity kernels of coda waves, in 2D and 3D, propagating in a medium with a realistic velocity model, were computed by numerical simulations following Yoshimoto [17]: assuming that the total seismic energy is subdivided in particles of unit energy that are randomly emitted from the source, the energy particle loses a fraction of its energy by anelastic absorption. When the energy particle encounters an elastic and isotropic scatter, it randomly changes its direction. After a number of collisions, it will arrive at the receiver at a given lapse time, measured from the origin time. Thus, the energy envelope is the sum of the energy particles arriving at the receiver from random directions at a given lapse time, sampling a portion of the Earth’s volume. Thus, it can be heuristically assumed that the more a volume element (in which the whole medium is divided) is crossed by energy particles, the more it weighs in the energy envelope of each seismogram. The space weighting function can be calculated in the whole propagation volume by counting the number of crosses and collisions of the energy particles in each volume element [11]. In the last few years, sensitivity kernels for scattering radiation have been used to image both volcanic structures [18,19,20] and tectonically active zones (e.g., [21,22]).

The first attempt to study the $Q_c$ attenuation in the Gargano Promontory (Northern Apulia, Southern Italy) (hereafter GP) has been performed using a dataset of microearthquakes recorded by the OTRIONS Local Seismic Network (OT), managed by UniBa, in the period from April 2013 to July 2014 [23]. The obtained 2D image of $Q_c$ revealed the presence of a well-resolved $Q_c$ anomaly in the southwestern sector of the GP, at $f_c>$ 5 Hz. The southwestern sector of the GP is characterized by the presence of fluid circulation in fractures at shallow crustal levels, as also revealed by 2D resistivity imaging [24] and by the heat flow map [25], which can cause both scattering and intrinsic attenuation effects on a local scale. The heterogeneity of the GP crust was recently highlighted by a local-scale seismotectonic study [26], thanks to a recently released dataset of microearthquakes recorded by both the National Seismic Network (IV), managed by the INGV, and the OT network in the period from August 2015 to August 2018 [27]. The authors [26] analyzed all possible focal mechanisms and found a great heterogeneity of solutions in the seismicity of the upper and intermediate crust; however, in the lower crust, the seismicity appears more concentrated in the NE sector, which is consistent with a thrust-faulting stress regime with maximum compressive axis ($S_{hmax}$) that is NW–SE oriented, and a shallow normal-faulting stress regime with a maximum horizontal compressive axis $(S_{hmax}),$ trending NW–SE [28].

The historical seismicity of northern Apulia reveals seismogenic structures that can generate destructive events. Among these, the 1627 earthquake (with an intensity X degree MCS and estimated magnitude of 6.7) and the 1646 earthquake (with an intensity IX–X degree MCS at Carpino and estimated magnitude of 6.2) [29]. Concerning the tectonic structures, many faults are recognized along the entire GP surface, the characteristics of which are poorly known. In Figure 1, we report the seismotectonic map of Italy ([30]) and the details on the GP sector, with the principal tectonic lineaments that could be responsible for the main historical earthquakes.

Considering these recent results [26], related to the strong lateral and vertical crustal heterogeneity of the GP, in this paper we performed a new and more robust study of $Q_c$ in order to address some still open questions: the first is the reason for the presence of lower crustal seismicity in the NE sector of the GP, with the absence of shallower seismicity; the second is the presence of very shallow seismicity in the SW sector of the GP, and its relationship with fluid circulation in the upper crust. This attenuation study benefits from the recently released dataset of GP microearthquakes [27,31] and from a 3D mapping approach based on sensitivity kernels [15]. The 3D tomography of $Q_c$ was obtained by combining the new dataset with the previous one [23], used for the 2D mapping of $Q_c$. With respect to the previous estimate of $Q_c$ in the GP [23], the used dataset is enlarged and consists of earthquakes with magnitudes ranging between 1 and 2.8, and the sampled area is extended to the NE, thanks to a change in the OT network geometry. All the details regarding the data and the network can be consulted in [27]. Moreover, 3D images of $Q_c$ for the GP, using a projection method based on weighting functions derived from sensitivity kernels [15] have been performed for the first time in this area. The present $Q_c$ tomography depicts the S-wave attenuation structure down to a 32-km depth, in three frequency bands centered at 3, 6 and 12 Hz, respectively. We will show how our present images provide evidence that the strongest attenuation contrasts are well correlated with the geotectonic setting of the GP.

2. Data and Methods

2.1. Dataset

The dataset consists of two hundred specifically located seismic events, extracted from the database of microearthquakes that has recently been released by [31] and described in [27]. We extracted data regarding the earthquakes recorded in the period from July 2015 to August 2018 to achieve better coverage of the GP crust, thanks to two OT stations that were added in the northern sector of the GP in June 2015. The OT network stations are equipped with short-period seismometers (Lennartz 3D-V seismometer (flat response above 1 Hz; 10 ms sampling rate)) [27]; the IV network stations are equipped with various types of seismic sensors, i.e., short period, broadband, and accelerometric instruments [32]. The selected dataset is composed of specifically located earthquakes ($RMS<0.2$ s; horizontal and vertical location errors of less than 1 km, magnitude $1\le M_L\le 2.8$) by using the 1D Gargano velocity model [33]. With respect to a previous study [23], we selected a dataset with a broader network geometry and a higher magnitude threshold, in order to validate the previous results obtained with very small earthquakes ($0.3\le M_L\le 1.8$) in a narrower region. Details on data acquisition and management are described in [27]. We selected 191 events occurring in the period from 1 July 2015 to 30 August 2018, recorded to a maximum of 21 seismic stations. Foci depths range between 1 and 35 km, with shallower events located in the SW, and deeper events located NE of the promontory, as already observed in previous seismological studies of this area [26,29,34,35] (Figure 2).

2.2. Coda Attenuation Method

$Q_c$ analysis as a function of frequency was performed by assuming the single isotropic scattering theory [1,38]. The theoretical decay rate of the coda energy envelope $A_{th}=A_{th}\left(r,f_c,t\right)$ is given by [39]:

$$\ln \frac{A_{th}r}{r_0\sqrt{K\left(a\right)}}=\ln A_0-\frac{\pi f_c}{Q_c\left(f_c\right)}{t}_L$$

(1)

where $r$ is the source to receiver distance (km), $r_0$ is the reference distance (km), $A_0$ is a source term, $f_c$ is the central frequency (measured in Hz), $K\left(a\right)=\frac{1}{a} \ln \frac{a+1}{a-1}$, $a=\frac{t_L}{t_s}$ with $t_s$ the S-wave travel time, and $t_L$ is the time (measured in s) elapsed from the earthquake origin time, with $t_L>t_s$. The observed energy envelope $A_{obs}=A_{obs}\left(r,f_c,t\right)$ can be computed from a seismogram $S\left(t\right)$ and its Hilbert’s transform $H\left(t\right)$:

$$A_{obs}\left(r,f_c,t\right)=\sqrt{H^2\left(t\right)+S^2\left(t\right)}.$$

(2)

The $Q_c^{-1}$ estimation for the investigated area was obtained by filtering the coda waves with a Butterworth two-poles filter in those frequency bands with a central frequency $f_c$ = 1, 2, 3, 4, 5, 6, 8, 10, 12, 14, and 16 Hz, with a bandwidth ranging inside the interval [$\frac{f_c}{\sqrt{2}}$, $f_c\sqrt{2}$]. Then, we least-squares fitted the envelope $A_{obs}$ versus the time using Equation (1), with $r_0$ = 1 km. With this procedure, $Q_c^{-1}$ can be estimated as a function of frequency $f$ at the assigned lapse time window $t_L$ [10]. In fact, even though the coda envelope decay rate is theoretically independent of both the source-receiver distance and of magnitude, it depends on the time $t_L$ elapsed from the event origin time. In addition, Aki and Chouet [1] assumed that scattering is a weak process and is not strong enough to generate secondary waves, once they encounter other scatters (Born approximation). Scattered waves are produced once seismic waves encounter heterogeneities, faults, cracks, or irregular topography [40].

Generally, the quality factor increases with frequency through a power-law relationship between $Q_c^{}$ and $f_c,$ due to Romanowicz and Mitchell [41]:

$$Q_c=Q_0{\left(\frac{f_c}{f_0}\right)}^{\alpha }$$

(3)

where $Q_0$ is a reference value of the quality factor at $f_0$ = 1 Hz, and $\alpha$ is the frequency dependence exponent, which is close to 1 and varies from one region to another region, on the basis of the heterogeneity of the medium [38]. Equation (3) indicates that the attenuation of seismic waves with the passage of time (i.e., the distance from the source) varies with frequency.

2.3. Coda Attenuation 3D Imaging

In order to obtain 3D coda attenuation images, we used the approach discussed by del Pezzo and Ibanez [11], based on the scattering kernels described by Pacheco and Snieder [15]. del Pezzo and Ibanez [11] computed the polynomial approximation of the analytically calculated sensitivity kernels, which noticeably reduces the calculation time costs. In particular, we used the sensitivity kernel as a weighting function, heuristically giving them a probabilistic meaning. The $Q_c$ weighting functions (del Pezzo and Ibáñez, 2020) are here denoted as:

$$w_{ij}\left(x_j, y_j, z_j,x_{si}, y_{si}, z_{si},x_{ri}, y_{ri}\right)$$

(4)

and represent the probability that, for the i-th source-receiver couple positioned, respectively, at {$x_{si}, y_{si}, z_{si}$} and {$x_{ri}, y_{ri}$}, the measured value of $Q_c$ effectively corresponds to its true value at the space point with the coordinates $\{x_j, y_j, z_j\}$. The subscript i spans from 1 to $N_{s-r}$ ($N_{s-r}$ is the number of source-receiver couples in the data set); j spans from 1 to $M_{pix}$ ($M_{pix}$ is the number of cells in which the space has been subdivided, here called “pixels”).

We indicated with $Q_{ci}$ the measure of $Q_c$ for the i-th source-receiver pair, and with $Q_{cj}$ the measure of $Q_c$ for the j-th pixel. In this heuristic scheme, $w_{ij}$ represents the un-normalized probability that the i-th measure, $Q_{ci}$, is associated with the j-th pixel. Averaging over the events and weighting for the kernel, at the single position j, the method yields:

$$Q_{cj}=\frac{{\sum }_i{Q}_{ci}{w}_{ij}}{{\sum }_i{w}_{ij}}$$

(5)

The measure set thus produces for any pixel, j, a set of point values of $Q_{ci}$ multiplied by their probability, $w_{ij}$. We consider their weighted average as the characteristic value of $Q_c$ at the specific point. In other words, we simply back-project the $Q_c$ measure by attributing the measure to a space point of given coordinates, as weighted for the weighting functions.

By applying the ordinary error propagation equation to Equation (5) and by assuming ${\sigma }_i$ the standard deviation of $Q_{ci}$, we obtain the standard deviation of the j-th pixel:

$${\sigma }_j=\frac{{\sum }_i{\left({\sigma }_i{w}_{ij}\right)}^2}{{\left({\sum }_i{w}_{ij}\right)}^2}$$

(6)

We used Equation (6) to estimate the error on $Q_{cj}$ in any j-th pixel. Since ${\sigma }_j$ represents the uncertainty of $Q_{cj}$ in the j-th pixel, and the diffusive waves are not related to seismic rays, we do not need to evaluate the ray crossing inside each j-th pixel. In fact, the weighting functions used are peaked at source and receiver. Consequently, the sole distribution of sources and receivers gives a sufficient idea of the areas where the coda radiation has the maximum probability of passing through. On the other hand, the coda tomography concept goes beyond the ray tracing assumption, taking into account the energy space distribution more than the wave front of the direct (ballistic) radiation, depicted by the ray trace space distribution. [20].

3. Results

3.1. Q_c Estimates

The collected seismograms were visually checked to prevent erroneous $Q_c$ estimates. The visual inspection of the envelopes was carried out on recordings filtered with a Butterworth two-pole filter with a central frequency $f_c$ = 8 Hz (bandwidth ranging between 5.7 Hz and 11.3 Hz). We removed all the recordings that present abrupt variations in the amplitude of the envelope, i.e., bumps of energy due to concomitant seismic events or other environmental or artificial sources. We also removed the recordings with sudden changes in the envelope slope in the portion of the seismogram immediately following the S-wave arrival [42]. The suitable recordings were processed with the addition of two time-markers (Figure 3):

• An initial marker, named $T3$, placed after the S-wave arrival, when the envelope shows a decreasing trend;
• A terminal marker, named $T4$, placed before a bump, or an abrupt change of the slope, or when the oscillations cannot be distinguished from the seismic noise.

The selected envelopes were linearly fitted inside the time interval [$T3,T4$] in a logarithmic diagram for different $t_L$ = 10, 15, 20, 25, 30, 35, 40 s and $f_c$ = 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 16 Hz, reaching a maximum of 70 measures of $Q_c^{comp}$ for each seismogram component. We accepted only the components for which the condition $T3<t_L<T4$ was satisfied. $Q_c^{comp}$ was computed from the slope $m$ of the linear regression in Equation (1): $Q_c^{comp}=-\frac{\pi f_c}{m}$. Any coda amplitude decay with a positive $m$ was rejected. In Figure 4, the number of $Q_c^{comp}$ estimates was plotted. It can be inferred that, for $t_L$ = 20, 25, and 30 s, the maximum number of $Q_c^{comp}$ estimates is reached at all frequencies, while for $t_L$ = 10 s the minimum number of $Q_c^{comp}$ estimates is available; all series show stability in the number of $Q_c^{comp}$ estimates starting at $f_c$ = 3 Hz, indicating that for $f_c=2$ Hz, the results are less robust.

All the analyzed envelopes, cut in the time window between the markers $T3$ and $T4$, are available in for the reproducibility of our results.

3.2. Q_c Averaged

To improve the quality of the $Q_c^{comp}$ dataset, we removed all those $Q_c^{comp}$ estimates with off-scaled values to avoid the invalidation of the average $Q_c$. The removed $Q_c^{comp}$ outliers represent 5% of the initial $Q_c^{comp}$ dataset. The error on each $Q_c^{comp}$ was computed from the error $∆m$ on the slope $m$ of the linear regression: $∆Q_c^{comp}=\left|\frac{\partial Q_c}{\partial m}\right|∆m$. $Q_c$ was then computed from a weighted average, by using the inverse square of $∆Q_c^{comp}$ as weights. The results are plotted in Figure 5, and it can be observed that $Q_c$ regularly increases at a fixed $f_c$ for an increasing $t_L,$ and at the same $t_L$ for increasing $f_c$, indicating a satisfying regularization of the data. The observed correlations between $Q_c$ vs. $f_c$ and $Q_c \mathrm{vs}. t_L$ are in total agreement with previous studies conducted in the same area [23] or in other tectonic areas as, for instance, in the central Apennines [39].

Taking the logarithm of Equation (3), we get:

$$\ln Q_c=\ln Q_0+\alpha \ln f_c.$$

(7)

Through the linear regression of $\ln Q_c$ versus $\ln f_c$, we obtained $Q_0$ from the intercept and $\alpha$ from the slope. We performed this analysis for all single-trace seismograms to obtain average $Q_0$ and $\alpha$ values for the different lapse times $t_L$. The results are reported in

Table 1. For each lapse time $t_L$, $Q_0$ and $\alpha$ estimates, with their errors and the radius of investigation in the Earth lithosphere, are reported.

${\mathit{t}}_{\mathit{L}}$	This Study	Radius of Investigation (km) 1
$10$ s	$Q_0$ = 11 ± 9; $\alpha$ = 0.96 ± 0.07	19.3
$15$ s	$Q_0$ = 24 ± 15; $\alpha$ = 0.80 ± 0.11	29.9
$20$ s	$Q_0$ = 31 ± 15; $\alpha$ = 0.88 ± 0.10	38.6
$25$ s	$Q_0$ = 34 ± 13; $\alpha$ = 0.97 ± 0.09	48.2
$30$ s	$Q_0$ = 40 ± 14; $\alpha$ = 1.00 ± 0.10	57.9
$35$ s	$Q_0$ = 49 ± 16; $\alpha$ = 1.00 ± 0.11	67.5
$40$ s	$Q_0$ = 63 ± 20; $\alpha$ = 0.95 ± 0.15	77.2

¹ Assuming the relations of Aki and Chouet [1].

, together with a comparison with other studies in the same area or in other worldwide areas. The increasing trend of $Q_0$ vs. $t_L$ in Table 1 reflects the size of the Earth’s volume, sampled by seismic waves. By adopting the simplifications of Aki and Chouet [1], it is possible to define an order of magnitude of the investigation depth. Estimating a mean value of $v_s$ = 3.86 km/s for the Gargano crust [33], we computed the radii of the spherical volumes (and therefore the maximum depths investigated) for all the considered lapse time windows (Table 1).

The average $Q_c$ for each station, named $Q_c^{st}$, is plotted vs. $f_c$ in Figure 6a–g. These plots exhibit the pattern of variation of the quality factor $Q_c^{st}$ with central frequencies $f_c$ at different lapse times $t_L$. The station $Q_c^{st}$ trend versus $f_c$ is irregular for $t_L$ = 10 s, for which the least number of $Q_c$ estimates is available; for $t_L$ > 20 s, more regular trends are observed.

In Figure 7, we plotted the $Q_c$ estimates with respect to other regions of Italy (Figure 7a) and other regions worldwide (Figure 7b). The same comparison is given in

Table 2. Comparison of $Q_0 \mathrm{and} \alpha$ of other regions with active tectonics in Italy and worldwide, for approximately the same $t_L$.

Region	${\mathit{Q}}_{\mathbf{0}}$	$\mathit{\alpha }$	${\mathit{t}}_{\mathit{L}}\mathbf{\left(}\mathbf{s}\mathbf{\right)}$
Gargano (Italy) [this study]	63	1	40
Gargano (Italy) [23]	50	0.7	40
Mts. Peloritani (Italy) [44]	50	0.7	40
Messina Strait (Italy) [45]	76	0.5	40
Southeastern Sicily (Italy) [46]	38	1.4	40
Central Appennines (Italy) [39]	88	0.6	40
Granada Zone (Southern Spain) [3]	86	0.9	40
Almeria Basin (Iberian Peninsula) [47]	63	0.9	40
Pannonian Basin (Croazia) [48]	108	0.7	40
Northern Morocco [49]	140	0.9	40
Racha Region (Georgia) [50]	79	1	60
Mainland Gujarat (India) [51]	112	0.9	40
Nepal Himalaya [52]	143	1	40
Northwestern Himalaya [53]	150	1	30
Sikkin Himalaya [54]	91	1	40
Taiwan [55]	93	0.8	30
Western Nagao, Japan [56]	112	0.7	30
Andaman Islands [57]	122	0.8	40
Parecis Basin (Brazil) [58]	98	1.2	40
Eastern Cuba [59]	64	0.8	20
North Tanzania [60]	146	1	40

for $Q_0$ and $\alpha$.

3.3. 3D Mapping

For the 3D mapping of $Q_c$, we followed the procedure described by [20] for the $Q_c$ tomography of Mount Etna. We considered for our imaging the results gathered at $t_L$ = 30 s. Three non-overlapping frequency bands, centered on $f_c,$ were selected: $f_c$ = 3 Hz ranging within [2.1, 4.2] Hz; $f_c$ = 6 Hz ranging inside [4.2, 8.5] Hz; $f_c$ = 12 Hz, ranging within [8.5, 17.0] Hz.

To make our images more reliable, we added to our results as hitherto described, as well as adding the results previously published by [23]. We then built a merged database of $Q_c^{comp}$ estimates and errors, $∆Q_c^{comp}$ (respectively, $Q_{ci}$ and ${\sigma }_i$ in Equations (5) and (6)). In this way, we collected a total number of source-receiver three-component estimates as follows: $N_{s-r}$ = 651 for each recording station, that is, 1953 $Q_{ci}$ values and ${\sigma }_i$ errors. The Earth medium was subdivided in $M_{pix}$ cubic pixels having side of 5 km. $Q_c$ images at each depth are representative of a layer of 5-km thickness. For the section maps, we used an 8-km depth step. Thus, we plotted $Q_c$ images at 0, 8, 16, 24, 32 km, each depth being representative of a 5-km j-th pixel. In Figure 8a–c, we plotted the $Q_c$ images for $f_c$ = 3 Hz, $f_c$ = 6 Hz and $f_c$ = 12 Hz, respectively. $Q_c$ tomographic images show very low errors since $\sigma$ is more than one order of magnitude less then $Q_c$ in the whole area at all the considered frequencies. Therefore, our results can be considered robust and significant.

4. Discussion

The quality factor $Q_c$ has been estimated in order to assess the effect of tectonic and seismic activity in the hazardous region of the GP.

Before drawing conclusions about the GP tectonics, the meaning of what $Q_c$ physically measures should be re-addressed. Within the frame of the single-scattering theory, $Q_c$ is considered representative of the total attenuation. The total attenuation of seismic waves is the sum of intrinsic $Q_i$ and scattering $Q_s$ attenuation: the first term accounts for the conversion of seismic energy into heat through anelastic absorption; the second term accounts for seismic energy redistribution through refraction, reflection and diffraction at random discontinuities distributed inside a homogeneous medium. The severe trade-off between the estimates of the intrinsic $Q_i$ and the scattering $Q_s$ is due to the mathematical formulation of the model itself [1] and it is impossible to deduce from the $Q_c$ estimates which quality factor ($Q_i$ or $Q_s$) prevails. Many studies indicate that, in a fairly homogeneous crust, $Q_i$ would prevail over $Q_s$ in the $Q_c$ value ([39,61,62] among others). This could be true on a regional scale; where an average Earth crust is sampled, heterogeneities appear as smoothed and the $Q_c$ distribution can be considered representative of the anelastic absorption. If $Q_c$ is analyzed with respect to frequency, it appears that at frequencies $f_c$ < 6 Hz, scattering is dominant, whereas, at frequencies $f_c$ > 6 Hz, intrinsic attenuation is dominant [9,63]. Monte Carlo simulations based on the law of energy conservation showed that for short lapse times, the single-scattering model is applicable, while for long lapse times, the diffusion model is more appropriate [64]. In particular, at long lapse times, when coda waves are in the diffusive regime, $Q_c$ asymptotically approaches the intrinsic $Q_i$ in a uniform half-space model [12].

By using different methods, the separation of $Q_i$ and $Q_s$ from $Q_c$ has been estimated both in Italy and worldwide, in many tectonic contexts: the Adaman Sea [65]; the Alborz region (Iran) [66]; the Almeria Basin, Southwestern Iberian Peninsula (for $f_c$ > 3 Hz) [47]; Spain and Turkey [67]; the Tokai region (Japan) [1]; Northeastern Italy [61]; the Messina Strait (Southern Italy) [45]. The results of all these studies suggest that the attenuation of coda waves is dominated by intrinsic attenuation. The same result was achieved for the Southern Apennines (Italy) [62], in a regional-scale study that included the GP.

In this work, with respect to a previous analysis [23], we used an improved dataset of microearthquakes that was one order of magnitude greater, recorded by both the OT and IV network stations, extending the seismic coverage of the area. The differences between the results of this study and those obtained in the previous one [23] are due to the different magnitude ranges considered. In fact, by increasing the event magnitude, the signal-to-noise ratio also increased, giving rise to more reliable results. Moreover, since the network geometry used for this study was larger than the previous one, the $Q_c$ values were averaged at different distances. Several observations from different regions worldwide [40] indicate an increase of the value of $Q_c$ when increasing the lapse time window $t_L$, due to the greater penetration of waves in the deeper Earth, where the seismic attenuation is less pronounced. Therefore, seven $t_L$ windows were taken, from 10 to 40 s, with a 5-s time step, to estimate the attenuation at different lapse times, i.e., at different depths. The penetration depth ranges between 19 and 77 km (Table 1). The frequency-dependence of $Q_c\left(f_c\right)$ in Equation (4) is interpreted as a tectonic parameter. Regions of high tectonic activity are characterized by low $Q_c$ when compared to stable regions.

Many measurements of $Q_0$ and $\alpha$, in different tectonic and geological settings, have indicated a trend of higher $Q_0$ (roughly 100 to 1000) and lower $\alpha$ (about 0.3 to 0.7) for stable continental regions, while a lower $Q_0$ (from 40 to 200) and a higher $\alpha$ (from 0.7 to 1.1) are found in tectonically active areas (e.g., [68,69,70]). In Table 2, we reported the results of the $Q_0$ and $\alpha$ measurements in several regions of active tectonics worldwide, for comparable values of $t_L$. In this study, $Q_0$ varied between 26 and 66 (Table 1), indicating that the upper lithosphere is seismically active. The frequency dependence exponent, $\alpha$, varies between 0.79 and 1 (Table 1); the lower value of $\alpha$ occurs for $t_L$ = 15 s, corresponding to a penetration depth of 30 km, in the lower crust, suggesting that the higher degree of heterogeneities corresponds to the depth where the seismicity of the GP is concentrated. $Q_c$, computed for each station, shows a low-quality factor and great variability at low $t_L$, while at high $t_L$, $Q_c$ estimates appear to be regular and higher (Figure 6). This result may indicate, as observed worldwide, that the upper crust is more heterogeneous than the lower crust, with the presence of high-frequency scattering. Stations of the OT network have very similar values and regular trends of $Q_c$ at all the lapse time windows $t_L$ (Figure 6). Regarding the IV stations, $Q_c$ estimates are always higher for stations that are inside the GP (SGRT, MSAG, APRC) than those outside the GP (AMUR, CAPA, MOCO, MRVN, SGTA, TREM) at all $t_L$, reflecting the heterogeneous characteristics of the crust in this area (Figure 6). It is worth noting that $Q_c\left(f_c\right)$ for the GP follows a trend very similar to that of the Central and Southern Apennines, Italy [39,62]. The Italian Apennines are the site of both historical and recent destructive earthquakes [71], so the similarity in $Q_c\left(f_c\right)$ with that of the GP may have evident and important implications regarding the seismic hazard. In fact, although for several decades, the GP has been affected by earthquakes of small or very small magnitude, it has historically been the site of very destructive earthquakes (Figure 1). Differences among $Q_c\left(f_c\right)$ estimated in this study and that estimated in other areas in Italy and worldwide (Figure 7a,b) can be ascribed to the different geological–structural settings derived from distinct geodynamic and paleogeographic domains.

Considering that $Q_c$ is representative of the intrinsic attenuation and that this similarity increases at high frequencies, as already discussed, we can interpret the tomographic images in Figure 8. All $Q_c$ images are coupled with the standard deviation $\sigma$, which indicates the high quality and robustness of the results. Seismic attenuation is homogeneous and high at all depths, at both low and medium frequencies (Figure 8a,b). A high $Q_c$ anomaly (marked with A in Figure 8a) is present in the southwestern GP sector at all depths but, since it corresponds to an increased $\sigma ,$ it is poorly resolved.

Some heterogeneities can be observed at high frequencies (Figure 8c). In the shallowest map sections (at 0 km, 8 km and 16 km) $Q_c$ shows a well-localized high $Q_c$ anomaly (labeled B in Figure 8c) and a north-eastward decreasing trend (toward label C in Figure 8c). The high $Q_c$ anomaly (zone labeled B, in Figure 8c) is observed at depth down to about 16 km in the northernmost part of the GP. This $Q_c$ anomaly is remarkable, since it corresponds with the area studied by Loddo et al. [72], where the author found the presence of a high-density and high-susceptibility body that was 10 km wide, at 4–7 km depth. The presence of this body agrees very well with the inferred low seismic attenuation in that portion of the crust. Loddo et al. [72] explained the presence of the deduced high-density/susceptibility body with the presence of an intrusive body within the sedimentary sequence, as outlined by the presence of minor outcrops of mafic intrusive rocks along the GP coastline. The observed northeastward-decreasing trend (toward the zone labeled C in Figure 8c) correlates well with the evidence of an anomalous absence of seismic activity in this sector at a depth of up to 20 km. In this frame, in the Northeastern sector of the GP, the low $Q_c$ in the upper/intermediate crust could be representative of ductile behavior, a speculation that deserves further investigation. We did not find any correlation with major surface fault lineaments of the area, as depicted in Figure 1.

In Figure 8c, at greater depths (24 km and 32 km), down to about 35 km, a high $Q_c$ anomaly (labeled D, in Figure 8c) appears in the easternmost sector of the GP. In this area, seismicity is deeper, as can be also observed in Figure 2. With respect to the previous 2D mapping of $Q_c$ [23], the easternmost part of the GP is covered by a seismic network, and it is now included in our results. It is very interesting that, in Figure 8c, in the two deepest map sections (24 km and 32 km), $Q_c$ shows a decreasing trend (toward label E in Figure 8c) toward the southwest: in this sector, at a depth greater than 18 km, seismicity is not present [26], the Adriatic Moho is shallower [73] and the crustal rheology has a ductile behavior [34]. All these independent observations agree well with the result of a lower $Q_c$ (zone labeled E, in Figure 8c). At the same depth, the easternmost sector of the GP is seismogenic [26] and in a fragile regime [34,73], and this agrees with the observed higher $Q_c$ (zone labeled D, in Figure 8c).

As discussed in a previous study [23], with increasing $f_c$, $Q_c$ shows the effect of the intrinsic and scattering attenuation of back-scattered body waves coming from deeper portions of the lithosphere [1]. At a higher frequency, in Figure 8c, it can be observed that, with increasing depth, $Q_c$ tends to decrease. This behavior can be ascribed to the effect of the intrinsic attenuation that increases with the sampled Earth volume. In confirmation of this behavior, as we plotted in Figure 9, the $Q_c$ trend with the source-to-receiver distance of $r$, with $f_c$ = 12 Hz and $t_L$ = 30 s, the linear decreasing trend of $Q_c$ ($r$) that confirms this hypothesis can be clearly observed.

In summary, from the $Q_c$ tomography at $f_c=$ 12 Hz (Figure 8c), it can be observed that in the upper/intermediate crust, $Q_c$ decreases toward the northeast (trend toward the C zone); in the lower crust, $Q_c$ increases toward the northeast (trend from the zone with label D to the zone labeled E). Therefore, the intermediate crust would act as a separation interface between these two trends of $Q_c$ in opposite directions. The attenuation maps confirm the high heterogeneity of the GP crust, as already evinced from the analysis of seismicity and the stress field [26,34].

5. Conclusions

Thanks to the improvement of the OT network, with the addition of two new stations in the northeastern sector of the GP, an enhanced dataset of microearthquakes was recently collected and released [27]. In this study, with this new dataset, the total $Q_c$ of S-waves was studied for the GP area. $Q_0$ values indicate that the Gargano area is tectonically active when compared with other regions as, for instance, the Central and Southern Apennines (Italy); the comparative results of $\alpha$ indicate that the GP crust is strongly heterogeneous at a penetration depth of about 30 km in the lower crust, and this result agrees with the observed seismicity, which concentrates at that depth, indicating a high degree of fracturing of the GP’s lower crust.

In this work, using attenuation kernels, we present the first 3D $Q_c$ maps of the Gargano (Southern Italy) area for different depths and frequencies. The geophysical knowledge of the area greatly benefits from this 3D approach, providing a marked improvement with respect to the averaged $Q_c$ estimates alone. The imaging was performed at different frequencies and depths. Results at low and medium frequency indicate a fairly homogeneous $Q_c$ at all depths, whereas, at high frequency, results indicate:

(1)

The presence of a high $Q_c$ body in the northern sector of the GP at a depth of up to 16 km, which would support the hypothesis of a high density/susceptibility body, as deduced by Loddo et al. [72] in the same sector and at the same depth (zone labeled B (Figure 8c));

(2)

The presence of a low $Q_c$ area, down to a depth of 16 km, suggesting a ductile regime in the upper/intermediate crust in the northeastern sector of the GP (zone labeled C (Figure 8c));

(3)

The presence of a low $Q_c$ area in the southwestern sector of the GP at depth down to 32 km that would support the presence of a ductile behavior in the lower crust [34], responsible for the absence of seismicity at that depths (zone labeled E (Figure 8c));

(4)

An increasing trend of $Q_c$ with depth (Figure 8c), confirming that, at high frequencies, $Q_c$ is dominated by the intrinsic attenuation of back-scattered body waves, as depicted in Figure 9.

According to these interpretations, the GP area reveals many interesting features that are promising for further investigation.