Structure of the crust and uppermost mantle beneath the Sicily Channel from ambient noise and earthquake tomography

In this work, we study the crust and the uppermost mantle structure beneath the Sicily Channel, by applying the ambient noise and earthquake tomography method. After computing cross-correlation of the continuous ambient noise signals and processing the earthquake data, we extracted 104 group velocity and 68 phase velocity dispersion curves corresponding to the fundamental mode of the Rayleigh waves. We computed the average velocity of those dispersion curves to obtain tomographic maps at periods ranging from 5 s to 40 s for the group velocities and from 10 s to 70 s for the phase velocities. We inverted group and phase speeds to get the shear-wave velocity structure from the surface down to 100 km depth with a lateral resolution of about 200 km. The resulted velocity models reveal a thin crust with thickness value of 15 km beneath the southern part of the Tyrrhenian basin and a thickness value of 20 km beneath Mount Etna. The obtained thickness values are well correlated with the reported extension of the Tyrrhenian lithosphere due to the past subduction and rollback of the Ionian slab beneath the Calabrian Arc. The crustal thickness increases and reaches values between 28 and 30 km beneath the Tunisian coasts and Sicily Channel. The S-wave models reveal also the presence of high velocity body beneath the island of Sicily. This finding can be interpreted as the presence of the Ionian slab subducting beneath the Calabrian Arc. Another high velocity body is observed beneath the southern part of the Tyrrhenian basin, it might be interpreted as the presence of fragments of the African continental lithosphere beneath the Tyrrhenian basin. of the study area with a spatial resolution length of 180 km. The Moho map reveals a shallow Moho with depth value of 15 km beneath the south of the Tyrrhenian basin, the depth value increases rapidly and reaches 25 km in the north-western part of the island of Sicily. This Moho topography likely reflects the oceanic-continental transition in the south-eastern part of the Tyrrhenian basin. In the Sicily Channel, we observe different Moho depth values in the northern part and southern part of the channel. The average Moho depth in the northern part of the channel is about 30 km, while it decreases and reaches a value of 25 km in the southern part. This is likely to be associated to the active NE-SW oriented extension reported by many authors in literature. The shear wave velocity model reveals also the presence of a high velocity body (Vs between 4.5 and 4.8 km/s) beneath the eastern part of the island of Sicily, this can be interpreted as the presence of the Ionian slab subducting beneath the Calabrian Arc. Other high velocity body is observed beneath the southern part of the Tyrrhenian basin and the northern and central part of the Sicily Channel, it might be interpreted as the presence of fragments of the African continental lithosphere. In conclusion, the tomography inversion of ambient noise and earthquake data allows us to provide useful informations about the seismic distribution and the Moho topography beneath the Sicily Channel and surroundings. The geometry of the crust and uppermost mantle, highlighted by this study, allowed us to have basic knowledge about the geodynamics occurring in the area that can be useful for further investigation.


Data and methods
The seismic ambient noise is defined as the constant vibrations of the Earth's surface at different seismic frequencies, even without earthquakes [Okada, 2003]. The seismic ambient noise tomography is considered as a passive tomographic method to image the Earth interior. Bensen et al. [2007], Shapiro et al. [2004Shapiro et al. [ , 2005 and Stehly et al. [2006Stehly et al. [ , 2008 showed that the cross-correlation of the permanent ground vibrations computed between a pair of receivers will result in a waveform that differs only by an amplitude factor from the Green function computed between the two receivers.
To study the velocity structure of the lithosphere and uppermost mantle beneath the Sicily Channel, we collected continuous recordings of ambient noise from 19 broad-band vertical-component stations within a time period of four years (01-01-2010 to 31-12-2013). Due to the limited number of available broadband Tunisian stations and to increase the number of rays crossing the area, we analysed 87 earthquake signals recorded at the Tunisian stations TAMR, TATN, and THTN, for a time period ranging from 02-04-2010 to 03-12-2016. The continuous recordings and earthquake signals were downloaded from different networks: Tunisia Broad Band network (TT), Italian seismic network (IV) (INGV Seismological Data Centre, [2006]) and Mediterranean very broadband seismographic network (MN) (MedNet Project Partner Institutions, [1990]), using the European Integrated Data Archive EIDA (http://www.orfeuseu.org/eida/eida.html). Map on the top right of Figure 1 shows the ray coverage of the study area.
To determine the surface wave dispersion of the Rayleigh waves fundamental mode from ambient noise data, we follow the processing data methodology proposed by Bensen et al. [2007]. For each continuous recording we (i) cut daily signals (ii) remove instrumental response that can obscure the ambient noise (iii) remove mean and trend (iv) apply a band-pass filter between 0.006 and 0.2 Hz (v) decimate the daily signal to 1 sample per second (vi) apply a time domain normalisation to reduce the effect of earthquakes on the cross-correlations, instrumental irregularities and nonstationary noise sources near each station (vii) apply a frequency domain by adding white noise to the data. We use the processed daily signals to compute the cross-correlation between each pair of stations and then we stack all the resulting daily cross-correlations to a single time-series to increase the Signal to Noise Ratio SNR. 68 cross-correlation signals with SNR greater or equal to 7 were retained. The retained signals are represented by red segments in the map on the top right of Figure 1. The processing and the computation of cross-correlations were done using the Seismic Analysis Code of Goldstein et al. [2003]. Figure 2b shows an example of positive and negative lags of the stacked cross-correlation between the two stations THTN and WDD depicted in Figure 2a, while Figure 2c shows the average signal obtained after summing the positive and negative lags. Figures 2d and 2e represent the group and phase velocity dispersion curves (red curves), respectively, that correspond to the Rayleigh waves extracted from the stacked cross-correlation signal. The dispersion curves are compared to two global models: Preliminary Reference Earth Model (PREM) developed by Dziewonski and Anderson [1981] (black dashed curve) and Global Dispersion Model (GDM52) (blue dashed curve) of Ekström [2011].
To increase the ray coverage, we used 87 earthquake signals occurred in Sicily from 02-04-2010 to 03-12-2016 and recorded on the vertical component of the Tunisian broadband stations: TAMR, TATN, and THTN. We selected 72 earthquakes with magnitudes greater than 4 and hypocentral distance greater than 100 km. The earthquake signals were processed by applying the same methodology from step ii) to v). The seismic travel times were automatically calculated using the automatic seismic travel time calculator TauP Toolkit [Crotwell et al., 1999].
To measure the dispersion characteristics of the fundamental mode of the Rayleigh waves from the ambient noise cross-correlations and earthquake signals, we apply a multiple filter analysis using the program Do_mft from Herman package [Herman, 2013]. We extracted 104 group velocity and 86 phase velocity dispersion curves. Figure 3 shows the example of group and phase velocity dispersion curves extracted from the cross-correlation of the ambient noise of 6 different station pairs.
To calculate the average speed of Rayleigh group and phase velocity dispersion curves at different periods on a 0.5° by 0.5° grid, we used the 2D tomography method of Ditmar and Yanovskaya [1987]. This method based on a generalisation to 2D of the 1D approach of Backus and Gilbert (1968), it has been widely used in surface waves tomography studies [Ritzwoller et al., 1998;Vuan et al., 2005;Guidarelly et al., 2011 andManu-Marfo et al., 2019]. It can also estimate the resolving power of the data by using a functional s(x,y) for different orientations of the coordinate system to determine the sizes of the averaging area along different directions, for each point (x,y) of the grid. The resolution can be approximated over the parameter L = (S min + S max )/2, where S min and S max are values of s(x,y), corresponding respectively to the largest and smallest axis values of an ellipse centred at that point [Yanovskaya et al., 1997]. The maximum period of each Rayleigh group velocity dispersion curve was fixed following  respectively. e) The red curve shows the Rayleigh-wave phase velocity dispersion picked using our data. The black and blue dashed curves are for the phase velocity dispersion obtained from the PREM and GDM52 models, respectively. The two models are well described in Ekström (2011). the criterion of Bensen et al. [2007], where the maximum period requires an inter-station distance Δ greater than three wavelengths ; (Δ > 3* = 3*c* ; where c is the Rayleigh group velocity). For each Rayleigh phase velocity dispersion curve, the maximum period was fixed following the study of Luo et al. [2015], where they demonstrated that Rayleighwave phase velocity dispersion measurements at short inter-station distances Δ as short as one wavelength are consistent and reliable as those at distance longer than three wavelengths.

Group velocity tomographic maps
Applying the method described above on our group velocity dispersion data, we get a spatial resolution of 180 km for periods from 5 s to 30 s. Figure 4 shows the resulting Rayleigh group velocity tomographic maps. The group velocity maps at shorter periods (5 s and 10 s) represent the sensitivity of the surface waves to shallow structures.

Phase velocity tomographic maps
Using the same tomographic method on our Rayleigh phase velocity dispersion data, we get a spatial resolution of 160 km for period 10 s, 180 km for periods 20 and 30 s and 250 km for periods 40, 50, and 60 s. Figure 5 shows the resulting Rayleigh phase velocity tomographic maps. The maps at periods 10 and 20 s are comparable to the respective group velocity ones. At longer periods between 30 and 60 s, the phase velocities are mapping the uppermost mantle materials beneath the study area, they are sensitive to velocity changes from 50 to 100 km depth. The velocity contrast is well correlated with the crust thickness. Negative velocity anomalies represent regions with thicker crust, while positive anomalies correspond to regions with thinner crust. The results are well correlated with crustal thickness values reported by Argnani [1990] and Marone et al. [2003]. The maps depicted in Figure 6a-f and 6g-l show the lateral resolution in km of the Rayleigh-wave group velocity and phase velocity, respectively, at each indicated period.

Shear wave velocity model
Analysing and plotting group and phase velocity derivatives with respect to elastic parameters versus depth, we found that group velocity periods are sensitive to S-wave velocity structure in the depth range approximately between 1 and 60 km, while the phase velocity periods are sensitive to S-wave velocity structure depth down to 100 km (see Figure 7).
To determine the shear wave velocity structure of the study area, we inverted simultaneously group and phase velocities, computed on each cell of the 0.5° by 0.5° grid used for the tomographic study. The method used for the inversion is the iterative linearized damped least-squares method proposed by Herman (2013 (Ludwig et al., 1970). To avoid the dependence on the initial models, we invert for each node of our grid, 20 different starting models using 5 different damping factors, letting the program repeat the inversion for 30 iterations.

Radia Kherchouche et al.
6 The starting models were parametrised as a stack of 1, 1, 24, 5 and 2 layers of 1, 2, 2.5, 5 and 10 km thick respectively and overlying a halfspace, in which the shear wave velocity was allowed to vary. The damping factor is the parameter used to find the best trade-off between the dispersion data misfit and the shear wave velocity model roughness. To choose the proper damping factors we follow the L-curve method [Hansen, 1992;1998], we retained the values of 0.01, 0.05, 0.5, 1 and 2.5. (see supplementary Figure S4). Due to the non-uniqueness of the solution, we obtained for each node of our grid 100 slightly different shear wave velocity models, the final solution is calculated as an average of all resulted models.
Analysing the resulted 1D shear wave velocity models, we determine the Moho depth value for each cell of the 0.5° by 0.5° using a velocity contour of 4.0 km/s as proxy. We interpolated the Moho depth values of each cell to get 7 Velocity structure of the Sicily Channel the 2D map representing the Moho topography beneath our study area (see Figure 8). in this eastern part of the Tunisian margins and the southern part of Sicily, this is well correlated with the low velocity anomalies revealed by the Rayleigh group velocity tomographic maps at periods of 5 and 10 s, which is interpreted as the presence of sediments in those areas. Other important feature is represented in the maps at 40, 50 and 60 km depth, by the presence of high velocity body (Vs between 4.5 and 4.8 km/s) in the eastern part of the island of Sicily, the southern part of the Tyrrhenian basin and the northern and central part of the Sicily Channel.
The final models resulting from the inversion method are plotted on 8 vertical sections representing the shear wave velocity structure of the crust and the uppermost mantle beneath the Sicily Channel. Figure 10 shows the obtained vertical cross-sections. Sections A, B, C, and D are crossing the Sicily Channel from northeast to southwest.
Sections E, F, and G are crossing the study area from the northwest to southeast. The cross-section H is along the 37° parallel (see the map in the bottom left of Figure 10). We plot also the seismic events from the catalog published

Results Discussion
In this study, we used the ambient noise tomography method to investigate the crust and uppermost mantle beneath the Sicilian channel. The tomographic maps at shorter periods (5 to 20s) reveal a negative anomalies in the southern part of the island of Sicily, the northern part of the Sicily Channel and the Tunisian margins, which is correlated with the presence of sediments at shallower depths of 5 and 10 km in the southern part of Sicily and in the eastern part of the Tunisian margins. These results are compatible with the Rayleigh-wave group and phase velocity tomographic maps, at the same periods, proposed by Manu-Marfo et al. [2019]. The presence of thick sediments between 5 and 15km is reported by Barberi et al. [2004]. and D'Agostino and Selvaggi [2004] from geological and geodetic data.
To study the resulted S-wave velocity models, we plotted 8 vertical sections mapping the crust and uppermost mantle beneath the study area. The vertical sections presented in Figure 10 with the plotted seismic events revealed very interesting features that help us to understand well the geodynamics of the Sicily Channel and its surrounding. to the Moho depths reported by Tesauro et al. [2010] from merging the most robust and recent Moho depth maps existing in the European region, and by Marone et al. [2003Marone et al. [ , 2004 from the joint inversion of local, regional and teleseismic data. The seismic activity observed in the Sicily Channel (see cross sections A, B, C, D, G and H of Figure   10) might be related to the rifting process affecting the channel and described by different authors [e.g. D'Agostino and Selvaggi, 2004;Serpelloni et al., 2007;Civile et al., 2008]. Other feature was revealed by cross sections C, D and G, were the presence of high mantle velocity values between 3.7 and 4.0 km/s is observed at depth between 15 and 20 km, beneath the southern part of the Tyrrhenian basin, and at depth between 30 and 40 km beneath the central part of the Sicily Channel. Similar high mantle velocity values are reported by Calò et al. [2013] using tomographic inversion of regional body wave phases. This seismic velocity distribution mimics well the oceanic-continental transition in the Tyrrhenian-Sicily domain. Other important feature is represented in sections B and C (Figure 10) by the presence of high velocity body (Vs between 4.5 and 4.8 km/s). This finding can be interpreted as the presence of the Ionian slab subducting beneath the Calabrian Arc. Another high velocity body is observed in sections F and G ( Figure 10) with similar velocities, this feature might be interpreted as the presence of fragments of the African continental lithosphere beneath the Tyrrhenian basin.

Conclusion
Submerged by sea water, the Sicily Channel remains one of the non-well studied areas in the central Mediterranean region, due to the limited station coverage on the North African side and the lack of OBS stations on the submerged part. In our knowledge, this study is the first attempt to determine the shear wave velocity structure beneath the Sicily Channel and surrounding, with spatial resolution length of 180 km for structures down to 50 km and with a resolution length of 250 km for structures down to 100 km. In this work we invert simultaneously the Rayleigh-wave group and phase velocities extracted from ambient noise and earthquake data, to determine the