Shear-wave velocity structure at Mt . Etna from inversion of Rayleigh-wave dispersion patterns ( 2 s < T < 20 s )

In the present study, we investigated the dispersion characteristics of medium-to-long period Rayleigh waves (2 s < T < 20 s) using both singlestation techniques (multiple-filter analysis, and phase-match filter) and multichannel techniques (horizontal slowness [p] and angular frequency [~] stack, and cross-correlation) to determine the velocity structure for the Mt. Etna volcano. We applied these techniques to a dataset of teleseisms, as regional and local earthquakes recorded by two broad-band seismic arrays installed at Mt. Etna in 2002 and 2005, during two seismic surveys organized by the Istituto Nazionale di Geofisica e Vulcanologia (INGV), sezione di Napoli. The dispersion curves obtained showed phase velocities ranging from 1.5 km/s to 4.0 km/s in the frequency band 0.05 Hz to 0.45 Hz. We inverted the average phase velocity dispersion curves using a non-linear approach, to obtain a set of shear-wave velocity models with maximum resolution depths of 25 km to 30 km. Moreover, the presence of lateral velocity contrasts was checked by dividing the whole array into seven triangular sub-arrays and inverting the dispersion curves relative to each triangle.


Introduction
High-resolution seismic images of crustal layers are extremely important for the retrieving of physical and mechanical properties of the materials.Knowledge of the crustal structure is of great relevance especially in a volcanic context, to detect magma intrusion and to improve the accuracy of earthquake location and source-mechanism definitions.The usual way to obtain an overall picture of the shear-wave velocities is through surface-wave dispersion analysis [Nur and Simmons 1969, Park et al. 1999, Xia et al. 1999, Dal Moro et al. 2008].The inversion of dispersion data for long-period surface waves has been extensively used to constrain the upper mantle shear-wave structure [Midzi 2001, Sabra et al. 2005, Shapiro et al. 2005, Cho et al. 2007].Moreover, many studies have obtained shear-wave velocity profiles from inversion of short-period (T < 3s) and mid-tolong-period (3 s < T < 20 s) Rayleigh-wave dispersion curves at different volcanoes, such as Popocatepetl [De Barros et al. 2008], Kilauea [Saccorotti et al. 2003], Masaya [Metaxian et al. 1997], Stromboli [Chouet et al. 1998, Petrosino et al. 1999, Petrosino et al. 2002], Vesuvius [Saccorotti et al. 2001] and Solfatara [Petrosino et al. 2006].To estimate the dispersion curves, these studies applied different techniques, including spatial autocorrelation analysis [Aki 1957], multiple-filter analysis (MFT) [Dziewonski et al. 1969, Dziewonski and Hales 1972, Herrmann 1987], and horizontal slowness and angular frequency (p-~) transform [Mc Mechan and Yedlin 1981, Herrmann 1987, Mokhtar at al. 1988, Russel 1988, Herrmann and Al-Eqabi 1991].
The aim of the present study was to improve our knowledge of the structure of Mt.Etna through surfacewave dispersion analysis.The most recent tomographic images of Mt.Etna are all based on the inversion of P-wave and S-wave arrival times from local earthquakes [Laigle et al. 2000, Chiarabba et al. 2000, Patanè et al. 2003, Patanè et al. 2006].However, the S-wave distributions that result from these inversions are poorly constrained, as S-wave picking in a volcanic context is difficult, due to the complex wave propagation in the typically highly heterogeneous and scattering medium that is found in volcanic environments.On the contrary, surface-wave dispersion curve analysis does not require phase picking, and allows for more robust determination of an averaged shear-wave velocity profile, and is thus less influenced by small heterogeneities.We have used mid-period and long-period Rayleigh waves from teleseismic (D [epicentral distance in angular degrees] > 13˚), regional (1˚< D < 13˚), and local (D < 1˚) earthquakes recorded by broad-band temporary arrays to determine the phase-velocity dispersion of Rayleigh waves in the frequency range from 0.05 Hz to 0.45 Hz.The inversion of the retrieved dispersion curves provides shear-wave velocity profiles down to a depth of 25 km to 30 km.Moreover, the existence of possible lateral variations was also investigated by dividing the network into sub-triangles and estimating the dispersion curves for each triangle.However, the results from this last analysis did not provide any clear evidence of lateral velocity variations within the region covered by our instruments.

The geophysical background of Mt. Etna
Mt. Etna is a 3,223-m-high stratovolcano that is located at the intersection of the main fault systems of eastern Sicily: (a) the NNW-SSE Aeolian Maltese fault, (b) the NE-SW Comiso-Etna-Messina system fault, and (c) the E-W Mt.Kumeta fault.Several geophysical and geochemical studies have been carried out to define the structural features of Mt.Etna volcano [Branca 2003, La Delfa et al. 2001, Catalano et al. 2004, Allard et al. 2006].The results of these studies excluded the presence of a huge magmatic chamber at a crustal level; on the contrary, some magma storage zones were located as embedded in the fracture systems.The main feature revealed by the most recent seismic tomography of Mt.Etna volcano is a high-Vp body that was interpreted as a solidified magma intrusion that extends from the south of the summit craters to the Valle del Bove, reaching a depth of 18 km [Villaseñor et al. 1998, Chiarabba et al. 2000, Laigle et al. 2000, Patanè et al. 2003, Patanè et al. 2006, Schiavone and Loddo 2007].The Vp values observed for this body are 5.6 km/s at 3 km in depth, and 6.5 km/s at 9 km in depth [Chiarabba et al. 2000].The high Vp body is embedded in several low velocity rings of sedimentary deposits that belong to the regional Apennine structure: flysch in the north-western and northern parts, and marine deposits in the southern part.The magma uprising is likely to be controlled by the structural weakness between the highvelocity body and the surrounding rocks, and probably occurs along the western and north-western border of the high-velocity body, where the main tectonic structures intersect.New and interesting data have been obtained from  the recent approach of four-dimensional (4D) tomography (3D velocity models, computed at different times).Using this technique, Patanè et al. [2006] identified a temporal variation in the Vp/Vs ratio during the flank eruption of 2002-2003, which they related to fluid injection in the upper layers that was possibly caused by a magma uprising.

Data analysis
The dataset we analyzed was recorded during two seismic surveys.The first lasted three months (October 31, 2002, to February 5, 2003), during the 2002-2003 eruption.Six Lennartz-MARSlite stations equipped with GURALP CMG-40T three-component broad-band seismometers (flat amplitude response curve in the 0.016-50 Hz frequency range) were installed during this period.The second survey lasted four months (from July to October, inclusive, 2005) and consisted of the installation of five Lennartz M24 stations equipped with Lennartz LE-3D/20s threecomponent broad-band seismometers, with a flat amplitude response curve from 0.05-50 Hz (Figure 1).Both of these seismic arrays were deployed and operated by the Istituto Nazionale di Geofisica e Vulcanologia, sezione di Napoli, using their mobile seismic network equipment.The minimum interstation distances were 2 km and 3 km for the 2002 and 2005 configurations, respectively.The array geometry in both cases was chosen as circular, to guarantee good azimuthal sampling of the waveforms.
For the surface-wave analysis, we used waveforms from teleseismic and local earthquakes with epicentral distances between 80 km and 12,000 km (Figures 2, 3 and 4).Moreover, the large distance range provided a wide frequency interval over which the dispersion curves can be estimated.We selected events with good signal-to-noise ratios (over a threshold of 10), measured in terms of the RMS of the signal SHEAR-WAVE VELOCITY STRUCTURE AT MT.ETNA and the noise, according to the relation: (RMSsignal-RMSnoise)/RMSnoise.We performed a polarization analysis to detect surface-wave arrivals, and then we applied singlestation and multichannel methods to the vertical components of these recordings, to obtain Rayleigh-wave dispersion patterns.The list of the events analyzed is reported in Table 1.

MFT and PMF analyses
To determine the local velocity structure from the inversion of the phase-velocity dispersion, array techniques that take into account only the traveling paths between the stations should be used.Single-station techniques are generally applied to determine group velocity dispersion curves that carry information on the whole of the medium between the source and the recording stations.
These are useful for preliminary analysis as they allow separation of the fundamental mode from the signal, suppressing the contribution of the higher modes.Therefore, before the application of the array techniques, we used MFT [Dziewonski 1969, Dziewonski and Hales 1972, Herrmann, 1987].This single-station method allows the determination of the group velocity dispersion through applying a series of narrow band-pass Gaussian filters centered at fixed frequencies that estimate the arrival time of the maxima of the surface-wave envelope at each central frequency.This procedure is based on the hypothesis that in a dispersed wave train, different components (at different frequencies) arrive at different times.The filter center frequencies were chosen in the 0.05 Hz to 0.45 Hz frequency range and the filter bandwidth was dependent on the epicentral distance.In this way, we extracted the preliminary group velocity dispersion curve for the fundamental mode from 17 earthquakes (Figure 5).
We also tried to apply MFT to extract the higher modes from the signals; however, in this case, the dispersion patterns were unclear for most of the seismograms, and therefore no useful information was retrieved from this MFT analysis of the higher modes.
We used the curve of the fundamental mode obtained for each event as the trial dispersion for the phase-match filter approach (PMF) [Turin 1960, Alexander and Lambert 1971, Herrin and Goforth 1977], which is an iterative procedure that is aimed at finding the filter that is phase-matched to a particular mode of a surface wave, thus allowing its extraction from the whole signal.In this way, for each event analyzed, we extracted the wave packet that corresponded to the fundamental mode, which was then analyzed using multi-channel techniques, including the p-~stack [Mc Mechan and Yedlin 1981, Herrmann 1987, Mokhtar et al. 1988, Russel 1988, Herrmann and Al-Eqabi 1991] and the cross-correlation methods, with the aim of computing the dispersion pattern of the phase velocity.

P-~stack analysis
The p-~technique provides an image of the phasevelocity dispersion curve in the slowness-frequency (p-~) domain by applying two linear transformations to the wavefield [McMechan and Yedlin 1981].We applied this method to the fundamental mode, previously recognized in the traces selected using the PMF method, and obtained the corresponding phase-velocity dispersion curves.These were finally averaged to retrieve the mean dispersion function in the frequency range from 0.05 Hz to 0.45 Hz (Figure 5).The average curve is indicative of the velocity structure of the whole area covered by the array.Possible bias in the phasevelocity dispersion curves can occur if the wave propagation deviates from the Great Circle Path (GCP), as the p-p rocedure does not take into account the deviation.Validation of the data obtained from the p-~stack in the 0.05 Hz to 0.12 Hz frequency range was carried out by applying an alternative approach, based on cross-correlation analysis.

Cross-correlation analysis
Under the hypothesis of a plane wave and that the diffraction effects out of the array should be removed [Pedersen et al. 2003, De Barros et al. 2008], we calculated the delay time of the surface-wave arrivals at the different stations from the time-lag at which the cross-correlation takes its maximum.We then inverted the estimated time delays using the L1 norm, and obtained the slowness and the propagation azimuth of the wave train that impinged upon the array.Through repeating this procedure in different frequency bands, we computed the phase-velocity dispersion curves for each earthquake.The errors affecting the phase velocity were calculated by propagation of the uncertainty associated with the delay time [Del Pezzo andGiudicepietro 2002, Menke 2002].
By averaging over all of the analyzed earthquakes, we computed the mean phase-velocity dispersion curve in the frequency range 0.05 Hz to 0.12 Hz (Figure 6).This procedure only gave good results for the seven teleseismic/regional events, as for the local earthquakes the array resolution was not enough to obtain a clear dispersion pattern.The wavelength for local earthquakes was about 3 km, whereas the array aperture (in terms of which the array resolution is defined) was about 15 km.In this case, there was a lack of coherency between the same signal recorded at different stations.As a consequence, comparisons between the results from the p-~and cross-correlation techniques were only possible in the 0.05 Hz to 0.12 Hz frequency band, where the dispersion patterns inferred from the two methods were similar, although the cross-correlation method provided smaller errors.On the contrary, in the 0.12 Hz to 0.45 Hz frequency range, reliable dispersion was obtained only from the p-~stack.

Inversion of dispersion curves and shear-wave velocity model
We inverted the average phase-velocity dispersion curves obtained from both the p-~stack and the crosscorrelation technique using a non-linear approach that was based on the neighborhood algorithm [Sambridge 1999, Wathelet et al. 2004, Wathelet 2008].This inversion procedure has the great advantage of being independent from the starting parametrization.To reduce the nonuniqueness in the inversion process, it is convenient to SHEAR-WAVE VELOCITY STRUCTURE AT MT.ETNA  fix the Vp values according to all of the previous information for the analyzed area.Therefore, we constrained the limits of the space parameter (Vp, Vs, and depth) based on the tomographic maps of Chiarabba et al. [2000], Patanè et al. [2003] and Patanè et al. [2006].The range of the starting Swave velocities was determined from the Vp values, allowing the Poisson ratio to vary in the 0.2-0.5 interval.We started the inversion procedure with the simplest space parametrization that accounted for the lowest number of layers, and then we increased it until we obtained the best fit between the observed and the predicted data [Scherbaum et al. 2003].This is the usual procedure used to solve the inverse problem with the lowest degree of freedom.The maximum and minimum solvable depths depend on the period range of the phase-velocity dispersion curve used in the inversion procedure [Dahlen and Tromp 1998].The inversion of both of the phase-velocity dispersion curves yielded similar velocity models with a maximum resolution depth of 25 km to 30 km (Figures 7 and 8).
Overall, the Vs models with misfit lower than 0.06 corresponded to theoretical dispersion curves that fell inside the error bars associated with the experimental dispersion.The Vs models were associated with lower misfit ranges, from 1.8 km/s to 2.5 km/s at shallow depths (1-3 km) and Figure 8 (bottom).Fits of phase-velocity dispersion curves from the cross-correlation method (left) and the Vs models (right).The Vs models obtained are compared with the Vs profiles extracted from previous tomographic studies on Mt.Etna.Vs profile labels as for Figure 7. from 2.8 km/s to 3.5 km/s down to 23 km.At depths greater than 23 km, the Vs reached a maximum value of 5 km/s, where the ensemble of models enlarged.The spread of solutions at greater depths can be related to the increasing uncertainty of the phase-velocity estimates at lower frequencies, to the decrease in the sensitivity with depth, and to the intrinsic nonuniqueness of the inverse problem.
The Vs models obtained by inversion of the average phase-velocity dispersion curves obtained from the p-~stack and cross-correlation techniques were very similar.This is easy to understand, as the two curves were perfectly overlapping within the error bars.However, the wide frequency range of the phase-velocity dispersion curve from the p-~stack gave better resolution for the shallower layers (1 km, as minimum resolvable depth), while for the phasevelocity dispersion curve from cross-correlation method, layers shallower than 5 km were hardly resolvable.
We compared the Vs models that resulted from the inversion procedure with some reference models extracted from tomographic studies of Chiarabba et al. [2000], Patanè et al. [2003] and Patanè et al. [2006].The 1D Vp models were extracted at various locations of the summit crater area inside the reference maps, and the corresponding Vs profiles were obtained by fixing the Vp/Vs ratio as 1.73.From the study of Chiarabba et al. [2000], we chose two models (vs_3km_W_Chiarabba00 and vs_5km_E_Chiarabba00 in Figures 7 and 8) for which the gradient of the velocity at depth was greatest.From the tomographic maps of Patanè et al. [2006], we extracted five profiles: one centered on the crater area, and the other four moving along the four cardinal directions from the center of the crater area, to a distance of 5 km (vs_5km_E_Patanè06,vs_5km_W_Patanè06,vs_N_Patanè06,vs_5km_S_Patanè06,vs_0km_Patanè06 in Figures 7 and 8).Moreover, from the tomographic study of Patanè et al. [2003], we extracted the Vp profile (vs_average_Patanè03 in Figures 7 and 8) averaged over all of the array area.Finally, we compared the theoretical phasevelocity dispersion curves that corresponded to all of these models to that obtained by resolving the direct problem using the best-fit model obtained by the inversion (Figure 9).All of the theoretically computed dispersion functions lie inside the estimated uncertainties of the observed phasevelocity curves.
To detect lateral velocity anomalies, we divided the whole array area into triangular subarrays [Laske et al. 2007], and in each triangle we computed the phase-velocity dispersion curve by the cross-correlation method associated with the multitaper spectral analysis [Park et al. 1987, Thomson 1982], to overcome the spectral leakage due to the use of single taper filters on seismograms.Application of the multitaper was limited to the analysis on triangular subarrays as there was no evidence of strong spectral leakage using the whole array dataset.For the triangular configuration, the number of available earthquakes was low, as in many cases there were less then three working stations.At frequencies lower than 0.06 Hz and higher than 0.08 Hz, the phase-velocity dispersion curves oscillated strongly, probably due to the poor resolution of the array and/or to local diffraction phenomena [De Barros et al. 2008].We inverted these dispersion curves, although in this case it was not possible obtain good resolution for Vs models shallower than 7 km in depth and deeper than 16 km.All of the Vs models found in the different triangular sub-areas were included perfectly in the Vs ensemble obtained by inverting SHEAR-WAVE VELOCITY STRUCTURE AT MT.ETNA the dispersion curves computed with the whole array; therefore, there was no evidence of any strong lateral variation in the 7 km to 16 km depth range.

Conclusions
The inversion of Rayleigh-wave dispersion curves obtained with different techniques has provided very robust S-wave (Vs) velocity models for the summit area of Mt.Etna volcano, which are compatible with recent models derived from seismic tomography.These velocity profiles that are derived directly from the dispersion curves are very reliable because the Rayleigh-wave velocity depends strongly on the S-wave velocity.The data from the inversion of the average phase-velocity dispersion curves from all of the arrays were fully compatible with the images of low (Vp/Vs) under the summit craters and the central part of the volcano [Laigle et al. 2000].This is supported by an indirect comparison between the Vs models obtained here and those retrieved from the Vp and Vp/Vs maps that have already been published by Laigle et al. [2000], with the extraction of an average vertical profile relative to all of the investigated area.The low S-wave velocities at shallow depths are compatible with the presence of fractured material, while the increasing Vs up to a depth of 18 km is in agreement with the presence of the high velocity body under the central craters [Chiarabba et al 2000].
The use of triangular sub-arrays with a horizontal resolution of the order of the triangle dimension (~10 km) does not allow the solving of any lateral velocity anomaly in the 7 km to 16 km depth range; however, we cannot exclude the presence of a small scale heterogeneity to which the Rayleigh waves are insensitive, as can be expected in a volcanic context.Moreover, the presence of possible lateral variations at shallower (<7 km) and deeper (>16 km) depths needs to be better investigated, by increasing the database and improving the array design in terms of the sizes and interstation distances, to gain greater resolution.A further extension to a 2D model of the S-wave velocity profile obtained in this study would allow for improved knowledge of the interior of the volcano, with possible implications for the definition of the location and size of a magma storage zone at crustal depths.

Figure 1 .
Figure 1.Map of the arrays used for the dispersion analysis, as installed in 2002 (blue triangles) and 2005 (red triangles).

Figure 2 .
Figure 2. Vertical components of a teleseism that occurred in Japan on 16.08.2005,recorded by the mobile network of the Istituto Nazionale di Geofisica e Vulcanologia (INGV), sezione di Napoli.

Figure 3
Figure 3(top).Locations of the teleseisms used for the dispersion analysis.Figure4(bottom).Locations of regional and local earthquakes used for the dispersion analysis.

Figure 5 .
Figure 5. Average phase-velocity and group-velocity dispersion curves obtained by applying the p-~transform and MFT, respectively, as indicated.

Figure 6 .
Figure 6.Comparison between the average phase-velocity dispersion curves obtained from the cross-correlation technique (cross-corr) and the p-~stack (p-omega) in the 0.048-0.10Hz frequency range.The error bars associated with the velocity values are standard deviations from the means of the earthquakes analyzed.

Figure 7
Figure 7 (top).Fits of phase-velocity dispersion curves from the p-~stack analysis (left) and Vs models (right) derived from the inversion procedure.The Vs models obtained are compared with the Vs profiles extracted from previous tomographic studies on Mt.Etna.Vs profile labels: vs_distance from the crater axis_direction_author of tomographic map; see text for further details.Figure 8 (bottom).Fits of phase-velocity dispersion curves from the cross-correlation method (left) and the Vs models (right).The Vs models obtained are compared with the Vs profiles extracted from previous tomographic studies on Mt.Etna.Vs profile labels as for Figure 7.

Figure 9 .
Figure 9. Comparisons between computed phase-velocity dispersion curves and those obtained solving the direct problem with the reference models from Chiarabba et al. [2000], Patanè et al. [2003] and Patanè et al. [2006].

Table 1 .
Analyzed earthquakes according to successful application of dispersion analysis.