Shallow velocity model in the area of Pozzo Pitarrone, Mt. Etna, from single station, array methods and borehole data

Seismic noise recorded by a temporary array installed around Pozzo Pitarrone, NE flank of Mt. Etna, have been analysed with several techniques. Single station HVSR method and SPAC array method have been applied to stationary seismic noise to investigate the local shallow structure. The inversion of dispersion curves produced a shear wave velocity model of the area reliable down to depth of about 130 m. A comparison of such model with the stratigraphic information available for the investigated area shows a good qualitative agreement. Taking advantage of a borehole station installed at 130 m depth, we could estimate also the Pwave velocity by comparing the borehole recordings of local earthquakes with the same event recorded at surface. Further insight on the P-wave velocity in the upper 130 m layer comes from the surface reflected wave observable in some cases at the borehole station. From this analysis we obtained an average P-wave velocity of about 1.2 km/s, compatible with the shear wave velocity found from the analysis of seismic noise.


Introduction
In volcanic environment the shallowest few hundred meters of rock are usually characterized by strongly variable mechanical properties. Therefore the propagation of seismic signals through these shallow layers is strongly affected by lateral heterogeneity, attenuation, scattering, and interaction with the free surface [e.g. La Saccorotti et al. 2001;García-Yeguas et al. 2011]. As a consequence, tracing a seismic ray from the recording site back to the source is a complex matter, with obvious implications for the source location [e.g. Ibáñez et al. 2000;La Rocca et al. 2004;Palo et al. 2009;Carmona et al. 2010]. For this reason the knowledge of the shallow velocity structure may improve the location of shallow volcano-tectonic earthquakes and volcanic tremor, thus contributing to improve the monitoring of volcanic activity [e.g. Saccorotti et al. 2004;Di Lieto et al. 2007;Mostaccio et al. 2013]. In particular, even if low velocity layers have been commonly observed on the top of the edifice at many volcanoes [e.g. Petrosino et al. 2002;Luzón et al. 2011;Cauchie and Saccorotti 2013], they are not usually considered in the localization analysis, inducing location errors. For all location methods in seismology, the greatest difficulty lies in finding consistent results as concerns the depth of the source, as for instance shown by Di Grazia et al. [2006] and by Zuccarello et al. [2013]. In fact, the procedure applied to the shallow volcano-tectonic earthquakes and volcanic tremor location, often relies on several assumptions, the most critical of which concerns the homogeneity of the propagation medium and the absence of site amplification factors. A deeper knowledge of the velocity model, especially in the shallow parts of the edifice, could mitigate the errors associated to the source location, and, eventually, the errors in the source mechanism analysis of the low frequency events (i.e. LP, VLP, volcanic tremor), as suggested by Bean et al. [2008] and De Barros et al. [2013].
Several methods have been developed to obtain a reliable estimation of the shear wave velocity in the shallowest layers. The most of them are based on the analysis of ambient noise vibrations, therefore they can be applied everywhere any time. The most used single station method consists in the evaluation of horizon-tal-to-vertical spectral ratio (HVSR) [Nogoshi and Igarashi 1971;Nakamura 1989], from which the resonance frequency of the local soil can be easily estimated. Information about the local velocity structure can be inferred from the frequency and width of the main HVSR peak [Wathelet et al. 2004[Wathelet et al. , 2008. The HVSR technique allows for a reliable detection of the fundamental resonance frequency fo that depends by subsoil stratigraphy, but it is less reliable in the estimation of local amplification values [Lermo and Chávez-García 1994].
The use of dense arrays of seismic stations installed at short distances (tens to hundreds of meters) has revealed a powerful tool for the investigation of shallow velocity structure [e.g., Aki 1957;Asten 1976Asten , 2003Asten , 2006Tokimatsu 1997;Scherbaum et al. 2003]. In fact array methods permit the estimation of the phase velocity and propagation azimuth of the analyzed signals. Determining the phase velocity as a function of frequency corresponds to estimate the dispersion of surface waves that propagate in a medium characterized by velocity changing with depth. Since seismic noise is rich of surface waves, and their phase velocity is related to the shear wave velocity, the analysis of seismic noise with appropriate methods and the successive inversion of the dispersion curve allows for an estimation of the shear wave profile of the shallow section. Many authors have studied seismic noise recordings with the aim of determining the shear wave velocity profile at shallow depth [e.g., Horike 1985;Tokimatsu 1997;Petrosino et al. 1999;Satoh et al. 2001;Scherbaum et al. 2003;Saccorotti et al. 2004;Wathelet et al. 2004;Di Giulio et al. 2006;Cauchie and Saccorotti 2013].
Aim of this work is the estimation of shear wave velocity profile in the area of Pozzo Pitarrone, NE flank of Etna volcano. We achieved our aim through the analy-sis of the dispersion characteristics of ambient noise vibration recorded by a temporary seismic array. This experiment is part the project MEDSUV funded by the European Commission (7th Framework Program -FP7 ENV.2012.6.4-2; grant agreement No. 308665), together with active seismic tomography experiment at MT. Etna -TOMO-ETNA [see in this volume: Coltelli et al. 2016;Díaz-Moreno et al. 2016;Ibáñez et al. 2016aIbáñez et al. , 2016b. We used the spatial auto-correlation method (SPAC) introduced by Aki [1957Aki [ , 1965 and later generalised by Henstridge [1979] and Cho et al. [2004], to estimate dispersion curves of surface waves. As described by Asten [2003Asten [ , 2006, the SPAC method has an advantage in processing seismic ambient noise because it extracts the scalar wave velocity of seismic energy propagating across the array, regardless of its direction [Asten 2006]. The insensitivity to the propagation azimuth is an important advantage when the sources of seismic noise are unknown.

Structural setting of Pozzo Pitarrone area
Mt. Etna is the largest stratovolcano in Europe and one of the most active and powerful basaltic volcanoes in the world. It is located in eastern Sicily, in a complex geodynamic framework where major regional structural lineaments seem to have a role in controlling the volcanic activity [e.g. Patanè et al. 2011]. Mt. Etna exhibits persistent voluminous degassing from the summit craters and recurrent summit and flank eruptions.
The data used in this study were collected in the area of Pozzo Pitarrone, which is a well drilled about one century ago for water caption, and now used for scientific purposes. The Pitarrone borehole is located in the middle northeastern flank, about 6 km from the summit craters, and close to one of the main intrusion zones of Etna volcano, about 1 km far from the so called NE-rift, and Pernicana fault [Azzaro et al. 2012]. The NE-Rift consists of sub-parallel eruptive fissures striking from 42°E to 47°E and extending from the surface (1700 m a.s.l) down to 2500 m depth. It is delimited eastward by a 200 m-high fault scarp in Piano Provenzana area. This fault, together with the Pernicana and Fiumefreddo faults, is commonly indicated as discrete segments of a near continuous left-lateral shear zone that dissects this sector of Etna (Figure 1b).
The eruptive fissures of the NE-Rift consist of spatter ramparts, scoria cone and coalescent scoria cones related to the eruptive activity of the last 15 ka belonging to the Mongibello volcano [Branca et al. 2011]. These products and the lava flows associated rest on the Ellittico volcano succession formed by rheomorphic lava generated during plinian eruptions occurred about 15 ka and by lava flows coming from the northern flank of the Valle del Bove [Branca et al. 2011]. The Ellittico lava flows consist of a thick succession that crop out in the area of Piano Provenzana forming the substratum in which the eruptive fissure of the NE-Rift are developed.

Instruments and data
Two broad band seismic stations of the permanent monitoring network are installed at Pozzo Pitarrone, one at surface (EPIT) and one at the well bottom, at depth of about 130 meters (EPIP). A temporary seismic array was deployed in the area around Pozzo Pitarrone between May and October 2014 ( Figure 1). The SHALLOW VELOCITY MODEL AT ETNA VOLCANO array consisted of 8 three component broad-band stations arranged in two roughly concentric rings of radii 100 m and 200 m centered at EPIT station ( Figure 2), deployed in an area doesn't cross by tectonic structures. The four stations of the internal ring were deployed at the end of May 2014, while the remaining group of stations was installed between June and July, 2014. In Table  1 the working status of the array stations during the period under investigation is shown. The seismometers used for the array were Guralp CMG-40T 60 s (http:// www.guralp.com) and Nanometrics Trillium 120 s compact (http://www.nanometrics.ca), and signals were acquired by Nanometrics Taurus data loggers recording at 200 sps, with 24 bits dynamic range (Table 2). During the acquisition period, the volcanic activity was characterized by continuous gas emission and occasional small explosions at the summit craters. Volcanic tremor was nearly persistent and significantly affected the recorded data.

Array response analysis
Any frequency-slowness signal processing method is based on a common waveform model of the signal. This means that each station of the array records the same waveform, except for a phase delay corresponding to the propagation across the array. Each plane wave component is assumed to cross the array with constant phase velocity, and each array station should record the wave with the same amplitude. These assumptions require that the seismic wavefield must be appropriately sampled both in space and time. Results of array analysis are reliable only in a limited frequency range, and depend on the number of stations, the array extension and configuration. Moreover, the spatial aliasing must be negligible in order that the solution be unique. The "array response function" is a good theoretical tool to evaluate the properties of the array configuration and the frequency range in which results of analysis are more reliable. We compute the array response by using the Beam-Pattern relationship in a version modified after Capon [1969]: where N is the number of sensors, x their position (vector x, y, z), w is the angular frequency, and S the slowness vector. For a given slowness S and angular frequency w, this function depends only on the relative position of the array elements, and gives the array response to a monochromatic plane wave with vertical incidence. The function B(S,w) given by Equation (1) is characterized by a peak of unit amplitude (normalized) centered at (0, 0) and many secondary peaks at different positions in the (Sx, Sy) space. The position and height of secondary peaks depends strongly on the number of stations and their relative position. The spatial aliasing is negligible when the height of secondary peaks is much smaller than one, and that happens when the array is composed of many stations placed at irregular positions to sample different distances and directions. The array response corresponding to the configuration used at Pozzo Pitarrone, computed for many different frequencies, is shown in Figure 3. Results indicate a small contribution of spatial aliasing in the array response. Regarding the frequency range, the broad peak observed at frequency 0.5 Hz indicates that the array resolution is poor at so low frequency, as expected being the signal wavelength larger than the array extension. The upper frequency limit is more difficult to estimate in practice because it depends on the station distances and by site effects, which reduce the coherence of the seismic wave-field among the array stations. In this case, the results of array analysis revealed to be reliable up to frequency of 4 Hz.

Data analysis
We analyzed the dispersion characteristics of ambient noise through the SPAC method, then we performed the inversion of the dispersion curves. The assumption for this kind of analysis is that the ambient noise is generated from a uniform distribution of sources. Assuming a constant seismic response (without local amplification) at each station of the array, the profile thus obtained provides a velocity model with one-dimensional characteristics. For this purpose we made a preliminary study of the local seismic response using the HVSR method.
First we performed the spectral analysis of the signals recorded by the array, and, successively, we analysed the spectral ratios between horizontal and vertical components (HVSR analysis) to find the occur- (1) rence of a resonance peak. This peak value was used in the following analysis, the estimate of 1D velocity model, in which we obtained a first shallow velocity model of shear waves. Finally, we performed the HVSR inversion analysis encompassing the characteristics to the 1D velocity profile in the Pozzo Pitarrone area.

Spectral characteristics and HVSR analysis
We have performed the analysis on the signal spectral components by using one hour of seismic signal that have the following characteristics: (i) constant ambient conditions, (ii) absence of volcanic activity, (iii) low level of volcanic tremor amplitude. Signal spectra are calculated by FFT over 20.48 s long windows (4096 sampling), sliding along the signal with 50% of overlap. The estimated spectra are characterized by a broad peak in the frequency band between 1 and 4 Hz at any array stations (Figure 4). In this frequency range the recorded signals are likely dominated by volcanic tremor characterized by amplitude peaks around 2 Hz. We achieved the spectral ratio analysis through the following procedure: 1. The windows of analysis were selected in order to maintain the stationary part of the signal and exclude transients; 2. A tapering was applied by using 5% cosine function; 3. Spectrograms were calculated by estimating the Fast Fourier Transform; 4. Spectra were smoothed by applying a logarithmic function given by the Konno and Ohmachi [1998] algorithm, which accounts for the different number of points with bandwidth smoothing of 40 samples; 5. The function H( f )/V( f ) was calculated for each window by the RMS mean; 6. Finally, the average HVSR function was calculated with its relative standard deviation.
A peak at frequency of about 0.8 Hz is evident at all stations ( Figure 5). It can be associated to a resonance effect even though the amplitude of this peak shows low quality if we consider the criteria suggested by the SESAME Project: amplitude greater than 2 and amplitude value of H/V ratio at the frequency f 0 smaller than half the amplitude in the range between f 0 and a quarter of f 0 [SESAME European project 2005]. The same amplifications are present at each station of the array, therefore they identify a site response constant in the area of investigation. Other peaks are present in the frequency range between 1 Hz and 4 Hz, probably produced by sources of volcanic origin like the volcanic tremor [e.g. Saccorotti et al. 2004;Ibáñez et al. 2009], as evidenced by the analysis of individual spectral components of the signal (Figure 4). Other peaks at higher frequencies (>10 Hz) likely represent resonance effects due to the presence of surface coverage widely distributed in the array area with variable thickness.

1D velocity model analysis
We used the software GEOPSY (http://www.geo psy.org) to estimate the propagation velocity of S-waves deduced from the surface waves dispersion properties. In the case of stationary noise, the cross-correlation method of Aki (SPAC method; Aki [1957]; Chávez-García et al. [2005]) is based on the properties of spatial SHALLOW VELOCITY MODEL AT ETNA VOLCANO correlation of the noise recorded at the vertical components of a circular array of stations with equal class distance r from a central sensor. Following the method described below, the dispersion curve of Rayleigh waves is used to determine the velocity structure. Let us represent harmonic waves of frequency ~by the velocity waveforms u(0, 0, ~, t) and u (r, i, ~, t), observed at the center of the array C(0, 0) and at point X(r, i) of the array. The spatial autocorrelation function is defined as: where u(t) is the average velocity of the waveform in time domain. The spatial autocorrelation coefficient t is defined as the average of the autocorrelation function z in all directions over the circular array: where t and i are the polar coordinates. If the signal is filtered in a narrow frequency band around ~0, the average autocorrelation coefficient, that is reported to the various pairs of receivers placed in different azimuth directions, is found by integration of Equation (3): where J 0 is the zero-order Bessel function of first order and c(~0), called dispersion function, represents the phase velocity at frequency ~0. From the SPAC coefficients t(r, ~), the phase velocity is obtained for every frequency from the Bessel function argument of equation, and the velocity model can be inverted [Rosas et al. 2011]. Many windows of signal with the same length are selected for each pair of array stations. The signals are filtered in 0.5 Hz wide frequency bands, and the correlation coefficients are calculated as a function of frequency, for each distance class. Each possible interstation distance is mapped to a circle of equivalent radius, and the set of such circles is called the co-array. Figure 6 shows the different rings chosen, with all station pairs.
The dispersion function c(~) is determined by equating the values found for the angular frequencies and the respective classes of distance to the arguments of the Bessel function J 0 . Using the analytical dispersion function it is possible to derive a 1D velocity model of the local surface structure. The inversion procedure requires a starting velocity model of the area, consisting of a set of plane-parallel layers overlying a half-space, each characterized by four parameters: thickness, S wave velocity, Poisson's ratio, density (h, Vs, v, t). The data inversion, carried out through the DINVER code (http://www.geopsy.org/), was performed through a high number of iterations that have produced about 5100 models, as shown in Figure 7a. The black dotted curve corresponds to the model that gives the minimum value of the misfit function, equal to 0.6, defined by Wathelet et al. [2004Wathelet et al. [ , 2008: (2) (3) (4) (5) in which V 0i represents the velocity at i-th frequency in the experimental dispersion curve, V RI represents the phase velocity of the fundamental mode of Rayleigh waves at i-th frequency. This function provides an estimate of the fit goodness of the experimental data to the resulting model. Figure 7a shows the resulting Vs profiles obtained by the inversion. Shear wave velocity between 500 m/s and 700 m/s characterize the shallow layer down to a depth of about 130 m. An abrupt velocity increase is observed at depth of about 60 m, sug-gesting the transition between two different layers.
We observe a good qualitative correlation between the estimated velocity model and the stratigraphic section of the investigated area. A lithological stratigraphy of Pozzo Pitarrone area has been inferred from detailed lithological description of the boreholes of the NE sector of Mount Etna (see Figure 7b) taking into account the different eruptive phases of the volcano history [Branca and Ferrara 2013, and reference therein]. In particular the borehole closer to Pozzo Pitarrone is en-SHALLOW VELOCITY MODEL AT ETNA VOLCANO  tirely made up of the products of the Stratovolcano phase [Branca and Ferrara 2013]. The Mongibello volcano consists of lava flow for the first 72 meters of thickness, separated from the Ellittico volcano products by about 15 meters of debris-alluvial deposit made up of volcanic clasts. By comparison between the stratigraphic thicknesses and the Vs profile obtained by our analysis, we find that the abrupt velocity increase matches the depth of the thick debris-alluvial deposit, as shown in the borehole stratigraphic column. In addition, the Vs values found in our analysis are consistent with a previous Vs profile calculated in nearby areas (Pernicana fault) through an active seismic experiment [Punzo et al. 2011].

HVSR inversion analysis
A technique that gives a more robust estimate of the Vs ground profile (shear wave velocity Vs-thickness of surface layers) is the joint inversion of the HVSR and the Rayleigh wave dispersion curves [Herak 2008]. Furthermore, inverting all observed HVSR curves, we obtain an extended two-dimensional characteristic to the velocity profile. The algorithm for the inversion estimates a set of parameters which minimizes the misfit function: Here indices OBS and THE stand for observed and the theoretical HVSR, and W i are the weights defined by Depending on W i , larger weights (for E> 0) are given to data around the frequencies where the observed HVSR is large. The model space consists of six parameters for each layer (excluding the half-space).
Since there are only three independent parameters in each layer for each wave type (P and S), and the density is common in both cases, the number of independent parameters in each layer is five. Uncertainties and bias of inverted values can be reduced by fixing some parameters (e.g., P-wave velocity or density in some layers in case they are sufficiently well-known from geotechnical explorations or other investigation methods). The joint use of the dispersion curve inversion technique, estimated by spatial autocorrelation, and inversion of the ellipticity curves, significantly reduces the level of non-uniqueness of the result. Even in this case the inversion of the data was achieved by DINVER code calculation, and it is relative to the . Shear wave 1D velocity profiles obtained by HVSR curves inversion. The black lines represent the velocity value corresponding to the minimum misfit. (6) same time window of signal noise used for the SPAC method. Figure 8 shows the velocity shear wave models relative to each station. The minimum value of the misfit function was obtained by a stratigraphic model characterized by three superimposed layers with increasing S-waves velocity. The first layer has a thickness of at most 4 m and is visible at the stations CT03, CT04, NA02 and NA03. The second layer extends to a depth between 50 and 60 m and it is characterized by a maximum S wave velocity of about 500 m/s. The transition between the second and the third layer is represented by a sudden increase of velocity at around 60 m depth.

P-wave velocity estimated from the borehole data
An estimation of the average velocity in the shallowest layer is obtained by the time lag measured between the signals recorded at the well bottom (EPIP) and at surface (EPIT). This estimate requires signals propagating nearly vertically between the two seismometers. Impulsive high frequency P waves are necessary for this kind of observation. We found some events appropriate for such measurements, as the one shown in Figure  9. The time lag between the direct P waves between EPIP and EPIT is about 0.1 s, as confirmed by the measurement of the time lag at the borehole station between the direct wave and the wave reflected from the free surface. This value correspond to an average P-wave velocity of 1.24 km/s when considering the way path of about 130 m. Through an active seismic survey curried out on EPIT station during the period under investigation, we obtained a Vp/Vs ratio of about 2. Assuming this value, from the shear wave velocity found from the inversion of the dispersion curve ( Figure 7a) we obtain Vp = 0.84 km/s in the first 58 m, and Vp = 1.40 km/s in the second layer. The corresponding average velocity is Vp = 1.10 km/s. The value 1.24 km/s found from the time lag between direct and reflected wave is about 15% higher than the average Vp value. This is not surprising for two reasons: first, the propagation of direct and reflected waves is not exactly vertical, thus the effective two-way path may be longer than two times the well depth; second, the well depth corresponds with the depth limit of our estimated velocity model, therefore it is possible that at such depth the effective velocity is a little higher than our estimation.

Discussion and conclusion
The analysis of seismic noise is a cheap but efficient tool for the investigation of the shallow velocity structure. In this work we took advantage from the availability of array data to estimate the shallow veloc-ity model of the area around Pozzo Pitarrone, NE flank of Mt. Etna. We applied the SPAC method to compute the dispersion properties of surface waves, then we inverted the dispersion curve to obtain a 1D velocity profile. To better constrain the inversion we used the HVSR computed at each array station, which also give a lateral 2D extension to the final velocity model. The obtained results indicate that site effects in the investigate area are quite homogeneous among the array stations. Our final result is a velocity model reliable to a depth of about 130 m, characterized by two main layers, the first one of thickness 60 -65 m and S wave velocity of 400-450 m/s, followed by a second layer with higher S wave velocity, estimated in the range 600-700 m/s. These velocity values are in good agreement with the properties of volcanic products that constitutes the shallowest layers in the NE sector of Mt. Etna. Moreover, the discontinuity found in the velocity model at depth of 60 -65 m is in good agreement with the stratigraphy of the area, as obtained from several boreholes located around the area of our investigation [Branca and Ferrara 2013]. Finally, the average shear wave velocity is in good agreement with the P-wave velocity estimated from the borehole data.
Studying the velocity structure in volcanoes is important for a better location of shallow sources. The next step in this direction is an investigation at greater depth. In fact studies like that presented in our work produce very detailed models, but they are always limited to depth of at most 200 m, but usually much less. Although it is reasonable that the shallowest 100 m of any volcanoes are characterized by similar velocity, the variation of velocity at depth from 100 m to 1 km are often the less known. In fact velocity models of the entire volcano estimated through tomography experiments often have a lower resolution, and near the surface present uncertainty or even gaps due to the stations distribution. Filling this gap of velocity models in volcanic environment could be achieved by the analysis of low frequency seismic signals recorded by broad band ar-SHALLOW VELOCITY MODEL AT ETNA VOLCANO Figure 9. A local earthquake recorded at 130 m depth (EPIP) and at surface (EPIT). The P wave reflected by the free surface (pointed by the arrow) is clearly recognized in the signal recorded at depth. rays characterized by tens of stations and extension of some km. The low frequency microseisms produced by sea waves could have the propagation properties appropriate to extend a velocity profile to depth of 1 km by using the same approach described in this paper.