sTOMO-ETNA was an active seismic experiment designed to improve understanding of the structure beneath Mt. Etna (Ibáñez et al., 2016a; Ibáñez et al., 2016b). The experiment involved the use of a large number of human and material resources (including seismic stations, and oceanographic, hydrographic, and logistical support vessels from Italy, Spain, and Greece). For this reason, it was decided to maximize the efficiency of the work by extending the shot locations from Ionian Sea to the Tyrrhenian Sea to include the Aeolian Islands (Coltelli et al., 2016).

The experiment, divided into several phases (between June and December 2014), generated active seismic signals using different capacities and powers of the air-guns. Those with the highest capacity (shots carried out in July 2014) were performed using compressors with 5,500 cubic inches of capacity. The aim was to generate two types of signal, one for producing images based on the refraction of seismic waves (using the maximum energy), and the other (based on lower power air-gun shots) to obtain images of marine structures using multi-channel techniques based on seismic reflection. Only shots generated with the maximum energy presented sufficient quality (based on the signal-to-noise ratio) to be used in tomographic structural models (both velocity and attenuation). Castro-Melgar et al. (2021) found that the most energetic shots had an equivalent duration magnitude of up to 1.4, based on the magnitude scale of Havskov et al. (2003) for volcanic regions. To this limitation we add a high content of background noise, which is typical of seismic signals in populated areas or on small islands (due to oceanic noise); as such, not all of the signals generated in the experiment were useful for this study. For our attenuation model, we only used signals generated in the Tyrrhenian Sea. From the set of all available seismic stations, we selected 37, guaranteeing at least one seismic station on each of the islands plus a number located on the island of Sicily (Figure 1). Seismic stations deployed on the island of Sicily (13 stations) belong to the portable seismic network provided by the Geophysical Instrument Pool Potsdam (GIPP), Germany. They are DATA-CUBE3 recorders and triaxial PE-6/B 4.5 Hz or Mark L-4C-3D seismometers. Seismic stations in the Aeolian Islands (24 stations) belong to the INGV (Istituto Nazionale di Geofisica e Vulcanologia) permanent network operated by the Osservatorio Etneo (the INGV Etna Observatory). Seismic stations of this network are broadband three-component Nanometrics Trillium seismometers (see Ibáñez et al., 2016a; Ibáñez et al., 2016b, for additional details). This distribution is not ideal for high-definition analysis, but was the best available. Of all the air-gun shots made in the Tyrrhenian Sea, we initially selected 3,450 (Figure 1). The available shot-station pairs provided an initial number of 127,650 waveforms.

The TOMO-ETNA database has previously been used to determine the attenuation structure of the Mt. Etna volcano region. For these studies, two methodologies have been used: 1) the separation of the intrinsic attenuation (Qi⁻¹) and scattering (Qs⁻¹) contributions through the diffusion model (Wu, 1985) and its 2-D spatial representation (Ibáñez et al., 2020); and 2) a 3-D attenuation tomographic study (Castro-Melgar et al., 2021) using the coda normalization method (Aki, 1980). Therefore, we tested the resolution of both of these methods in this study. The coda normalization method was ruled out owing to the poor resolution of the results (based on checkerboard and isolated anomaly tests). However, since the theoretical basis for the representation of the separation of Qi⁻¹ and Qs⁻¹ assigns attenuation values to a broader spatial region (e.g., Del Pezzo et al., 2018; Del Pezzo and Ibáñez, 2020), the checkerboard test indicated that a significant part of the region to be studied could be interpreted in a reliable way using this approach.

The diffusion method to obtain Qi⁻¹ and Qs⁻¹ and its spatial representation are explained in detail by numerus authors (e.g., Wegler and Luhr, 2001; Del Pezzo, 2008) and have been widely used (e.g., Prudencio et al., 2013a for Tenerife Island; Prudencio et al., 2013b for Deception Island; Prudencio et al., 2015c for Stromboli volcano; Prudencio et al., 2017a for Asama volcano; Prudencio et al., 2017b for Usu volcano; Prudencio et al., 2018 for Long Valley). This method is based in the assumption that the seismogram energy envelope can be modelled by the so-called radiative transfer equation (also called the transport model; Sato et al., 2012). Wegler and Luhr (2001) demonstrated that when the propagation medium can be assumed as extremely heterogeneous, as in volcanic structures or the shallowest part of the crust, then the diffusion approximation of the transport model is valid. In this approximation, the seismogram energy envelope is developed as a function of lapse time and source–receiver distance in terms of intrinsic and scattering attenuation coefficients. In this approach the seismic energy density as a function of time and source–station distance is linearised with respect to the two attenuation parameters, making it possible to separate the contributions of scattering and intrinsic attenuation effects through simple linear inversion (Prudencio et al., 2013a). Here, we simply highlight that by using the diffusivity (d) and coefficient for intrinsic attenuation (b) to estimate the intrinsic and scattering attenuation Q factors (Q _i and Q _s , respectively) as follows: Q i = 2 πf b (1) Q s = 2 πfpd v 2 (2) where f is the frequency and v is the velocity of S waves. In addition, the total attenuation value (Q_t ⁻¹) is represented as: Q t − 1 = Q i − 1 + Q s − 1 (3)

As described in the next section, the determination of the coefficients b and d is carried out by fitting the energy of the entire seismogram, after the arrival of the P-wave, according to the equations representing the diffusion model.

The data analysis process was systematically structured to ensure the highest quality of data. This process is summarised in the following steps.

To plot the spatial distribution of attenuation parameters, we followed the approach of Del Pezzo et al. (2016), as also described in Del Pezzo and Ibáñez (2020). Ibáñez et al. (2020) provide a detailed description of this procedure for data associated with Mt. Etna volcano. In brief, a back-projection method is used to plot spatial variations of Q _i and Q _s using a space-weighting function. We numerically estimated the energy envelope using the Monte Carlo approach of Yoshimoto (2000), who assumed that energy particles propagate according to Fermat’s rule and change direction on meeting randomly distributed scatterers; this distribution depends on the attenuation parameters. At the receivers, the sum of their energies defines the energy envelope; as such, the weighting function can be estimated. We obtained an intrinsic attenuation weighting function from the path density, and a scattering attenuation weighting function from the collision density. For each source-station pair, the Q _i and Q _s distributions were obtained by applying the weighting function described in Del Pezzo et al. (2016). We divided the area into cells, and for each cell we assigned corresponding Q _i and Q _s values. The final spatial distribution corresponds to total, intrinsic, and scattering values obtained by averaging all values assigned to every cell. This method is dependent on the experimental configuration, and is not affected by data quality (Prudencio et al., 2017a,; Prudencio et al., 2017b).

Following Ibáñez et al. (2020), we performed a series of tests to determine the quality of the spatial representation of the Qi, Qs, and Qt values. In contrast to the ordinary checkerboard-like 2-D test, as described in Prudencio et al. (2013a) and Prudencio et al. (2013b), we used a spatial representation of the uncertainty associated with the estimates of Q values in each pixel, as described in detail in Del Pezzo and Ibáñez (2020). This method was the most suitable for the present case, in which we did not perform an inversion of parameters as in most seismic tomography studies, but a projection mapping of the estimated parameters in the space through the weighting functions. To complete the uncertainty representation, we also present the density of the spatial ray path coverage (Figure 4), where a higher density of rays indicates better represented areas. In Figure 5 we show the pattern of “‘Resolution” quantity (see Del Pezzo and Ibáñez, 2020) and the standard deviation associated with the single Q values in each pixel.

The area with the best ray coverage (Figure 4) presents the best accuracy in our results, and is associated with a region covering the southernmost islands, the portion of the Tyrrhenian Sea between them, and northern Sicily. The entire archipelago of the Aeolian Islands is well resolved; the area located to the northwest of the selected region has the poorest reliability, although the errors do not exceed 25%. As seen in successive figures, the spatial distribution of the attenuation anomaly implies variation of up to 120%; as such, although the values themselves may not be exact, interpretations based on deviations from the mean will be reliable. Therefore, we conclude that the region selected for study is well represented by our seismic attenuation maps.

Since our approach provides single Q _t , Q _s , and Q _i values for each shot-station pair, lapse time, and frequency band, it is possible to obtain average values for the region. It should be noted that in the study of seismic attenuation there is an indistinguishable, but sometimes confusing, use of the terms Q and Q⁻¹. Here, seismic attenuation refers directly to the term Q⁻¹, since the attenuation is directly proportional to Q⁻¹. The parameter Q represents the so-called quality factor; that is, it is the ratio of total energy to energy loss, which is inversely proportional to attenuation. When we refer to numerical seismic attenuation values, it is common to speak in terms of Q values, since Q takes positive values greater than 1. In general, highly attenuating structures are reflected by low Q values and vice versa. We obtained the average values for Q _i ⁻¹ , Q _s ⁻¹ , Q _t ⁻¹ , Q _i , Q _s , and Q _t for lapse times of 20 and 30 s (Tables 1, 2, respectively).

We found that scattering dominates over intrinsic attenuation for all frequencies. As such, seismograms reflect geological heterogeneities at shallow depth. Using local and regional earthquakes with a source-receiver distance of 5–70 km, Del Pezzo et al. (2019) analysed intrinsic and scattering attenuation for a region that included our study area using the Multiple Lapse Time Window Analysis (MLTWA) method (Hoshiba, 1991; Akinci et al., 1995). They also found that scattering dominates over intrinsic absorption in the Aeolian Islands region. However, while we observed this for all frequency bands, they found that at higher frequencies intrinsic attenuation was stronger than scattering. These differences reflect differences in the depths of seismic sources. In this study, seismograms were generated at the surface (air-gun shots) and provide data pertinent to the uppermost crust (i.e., the first few kilometres). In Del Pezzo et al. (2019), earthquakes were deeper (∼10 km) and could be recorded at distances of up to 100 km, and thus provided data to resolve deeper structures. These differences confirm that the Q factor increases with depth (Ibáñez et al., 1990; Badi et al., 2009).

Figure 6 compares Q _i and Q _s values obtained for different volcanic areas in studies performed by our research group. We compared results using similar time lapses; the shortest (15 s) was used for Deception and Stromboli islands and the largest (30 s) was used for the Long Valley case. For the remaining areas, a 20 s lapse time was used. We observe that Q _s is always significantly lower than Q _i , implying that scattering attenuation is stronger that intrinsic attenuation for all compared areas. Unsurprisingly, results derived using the MLTWA method (using earthquakes of different depths; Del Pezzo et al., 2019) show much higher Q (i.e., lower seismic attenuation) values (for both scattering and intrinsic attenuation). However, when comparing only those values obtained using the same method and data type, we generally observe the following results:

1) Q _s values for the Aeolian Islands exceed those of other volcanic regions. We suggest that this difference reflects subsurface structural characteristics. Although the region is obviously volcanic, the results are also influenced by the attenuation of oceanic crust, which has lower apparent scattering.

These results justify the need to obtain maps of the spatial distributions of both scattering and intrinsic attenuation in order to determine how volcanic structures affect these values.

Spatial representation of attenuation was performed by plotting the distribution of Q⁻¹ values (directly proportional to the attenuation) rather than by plotting the absolute values. The Q⁻¹ values were calculated as the difference in percentage from the average value, with a maximum of up to ±120%. We use a rainbow colour scale, where the average value of the seismic attenuation is represented in green, areas of low attenuation are represented by cold colours (tending towards blue), and areas of high attenuation are represented by warm colours (tending towards red). With the exception of some peripheral areas, the study region is well-resolved.

As previously mentioned, the mean attenuation by scattering dominates over the intrinsic attenuation. We investigated whether this trend is general for the entire region under study or if there are structural elements for which the opposite occurs. Figure 7 shows anomaly maps for intrinsic, scattering, and total attenuation (represented by Q⁻¹) for the 4 Hz frequency of the two lapse times analysed (20 and 30 s, as shown in Figures 7A,B, respectively). In these maps we have used the same colour scale for all values, for Q _T , Q _S , and Q _i . On this scale, the lowest Q values are represented in red and the highest in blue. In this way, areas with the greatest attenuation are dominated by red hues, and those with the lowest attenuation by blue hues.

In general, scattering attenuation dominates over intrinsic attenuation throughout the region. To better illustrate this observation, Figure 8 shows the Q _S /Q _i ratios for the 20 s time lapse and the frequencies of 4, 6, 8, and 12 Hz. As observed for the entire region, this ratio is less than 1, which confirms that attenuation by scattering clearly dominates over intrinsic attenuation. We must emphasize that around the islands of Stromboli, Lipari, and Vulcano, this scattering effect is an order of magnitude greater than that of areas located further to the west, where marine deposits predominate and where there is a certain similarity between the two attenuation mechanisms.

We made maps of the spatial distributions of intrinsic (Figure 9) and scattering attenuation (Figure 10), from which we identified a set of attenuation anomalies. We interpret as anomalies those zones that differ significantly from the mean value of the region (represented by green hues). From the intrinsic attenuation maps, we identified at least seven regions (labelled as I1 to I7) and from the scattering maps another four anomalous areas (labelled as S1 to S4).

Table 3 shows the average size of the main identified attenuation anomalies, both for intrinsic and scattering attenuation. This sizes vary between a minimum area of ∼400 km² to some of more than 2,400 km². As the total size of the studied region is 16,800 km², these represent proportionally small-scale attenuation anomalies.

At this time we do not have an explanation for all of these regions; however, for the others we are able to offer some interpretation.

In reference to the spatial distribution of the intrinsic attenuation (Figure 9), we offer the following observations:

(1) From Tables 1, 2, there is not a strong dependency on frequency. In areas of high intrinsic attenuation, the dependence on frequency is even lower.

As discussed, this high intrinsic attenuation likely reflects sedimentary deposits on the seabed, which in this region are enhanced by the additional presence of volcanoclastic deposits (e.g., Romagnoli et al., 2013). According to De Ritis et al. (2010), these areas are also characterised by negative density and magnetization contrasts owing to the thick sedimentary cover.

In reference to the spatial distribution of the scattering attenuation (Figure 10), we make the following observations:

(1) Attenuation by scattering is more dependent on frequency than intrinsic attenuation. In this way, an evident shift from high attenuation anomalies at low frequencies to low attenuation anomalies at high frequencies is observed.

It is possible to explain these two observations based on the volcanic nature of the region; that is, on the presence of volcanic islands, where the dependence on frequency is obvious. As discussed, the islands are small in size, and scattering depends on the relationship between the dimensions of the heterogeneity and the wavelength of the incident wave. The fact that scattering attenuation tends to disappear at high frequencies (i.e., at shorter wavelengths) shows that it is the volcanic structures of the islands as a whole that most influence this type of seismic attenuation. For this reason, scattering attenuation is more relevant than intrinsic attenuation. Finally, the two regions marked as S1 and S2 are generally associated with volcanic structures and if these results are compared to the velocity tomography of Díaz-Moreno et al. (2018) there is a high concordance between high seismic attenuation by scattering and low velocity of propagation. In other words, structures of magmatic origin produce greater scattering of seismic waves. It is interesting to highlight that these two areas of high attenuation by scattering were also identified in a pioneering study of attenuation carried out by Del Pezzo et al. (1979).

On the other hand, at the time lapse of 20 s (Figure 10) there is a small region located to the west (marked as S4) where there is a low scattering effect and medium or low intrinsic attenuation. As this region is outside of the volcanic zones and is only observable for a 20 s time lapse (i.e., for the most superficial areas) we interpret this as attenuation due to the presence of marine sediments in the oceanic crust. This area belongs to the Cefalù Basin, a structural depression filled by a sedimentary cover formed by Tortonian to recent deposits (Pepe et al., 2004). The thickness of the Cefalù Basin cover is variable and can reach 1.5–2.0 km (Pepe et al., 2004; Milia et al., 2018), and gives rise to negative gravimetric and magnetic anomalies (De Ritis et al., 2010). According to Milia et al. (2018), the sedimentary cover of the Cefalù basin is made of deep-water turbidite deposits.

In summary, our results show that this area is not a large volcanic structural complex such as that of nearby Mt. Etna volcano, but rather it contains multiple dimensionally smaller volcanic structures.

It is also worth noting that the S1 feature spatially coincides with the northernmost part of the so-called Aeolian–Tindari–Letojanni fault system (e.g., Ghisetti, 1979; Lanzafame and Bousquet, 1997; Palano et al., 2012; Barreca et al., 2014). This system extends from the central sector of the Aeolian Archipelago to the Ionian coast and is characterised by right-lateral kinematics (Alparone et al., 1998; Ventura et al., 1999; De Astis et al., 2003; Ventura, 2013). Moreover, the unusual NNW–SSE orientation of the islands composing the central sector of the Aeolian archipelago (Salina, Lipari, and Vulcano), in contrast to the general geometry of the arc, is considered to be the expression of such a major regional tectonic structure, whose interpretation is still open to debate (Barreca et al., 2014 and references therein). It is not surprising that a rock volume affected by a fault system is characterised by high seismic attenuation by scattering. Based on a 3-D distribution of scatterers in central California, Nishigami (2000) found that the scatterer distribution was roughly correlated with the surface trace of the San Andreas fault system. High attenuation and high-scattering anomalies were also reported by De Siena et al. (2014) for fractured volumes of the Mt. St. Helens Seismic Zone.

In the present work we have demonstrated that kernel-based, separated intrinsic- and scattering-attenuation imaging allows the detection of local lateral contrasts in attenuation even for small lateral contrasts. This technique offers new opportunities study complex regions, including volcanic regions dominated by small islands or several isolated volcanic complexes, in order to better constrain their structures and improve understanding of their dynamics.