Wavefront distortion across Mt . Vesuvius observed by a seismic array

Local and regional earthquakes recorded at Mt. Vesuvius in 1997 and 1998 by a temporary seismic array have been analyzed with array methods to estimate slowness and backazimuth of correlated phases. The backazimuth observed for the direct P-wave has been compared with the expected direction, and results show surprising large differences in most of the cases. Since the observed difference is not systematic but depends on the epicenter backazimuth, we ascribe it to a wavefront distortion produced by the propagation through a non-homogenous velocity structure. To check this hypothesis we computed synthetic wavefront propagation in two different velocity models by applying a finite-difference method. The first model was homogeneous, while the second one was heterogeneous, with 2D symmetry in the horizontal plane given by a positive velocity anomaly centered at the volcano crater. We observe that the heterogeneous model gives results in good qualitative agreement with the observations. Such model is supported by the results of tomographic studies carried out for Mt. Vesuvius.


Introduction
Heterogeneity in the Earth's crust span a wide range of size and impedance contrast.Beside the generation of scattered waves, they are responsible of ray path and wavefront distortion in a broad range of wavelength.These phenomena are observed in a variety of environments and in a broad frequency range, from low frequency teleseismic phases to high frequency signals from local sources.
An useful tool to detect and analyze the distortion of seismic wavefront is given by the array techniques that provide detailed information about the propagation direction and velocity of coherent signals thanks to the spatio-temporal sampling of the seismic wavefield.Many studies focused their attention on tectonic areas [Arlitt et al. 1999, Shulte-Pelkum et al. 2003] by using array deployments, and related the slowness and backazimuth anomalies to crustal discontinuities beneath the recording sites.Small scale heterogeneous structures are very common in volcanic regions and the effects of the propagation throughout such media are described in many papers.Garcia Yeguas et al. [2011] analyzed slowness and backazimuth anomalies at Deception Island volcano (Antarctica) and interpreted the results as produced by the propagation through a shallow magma chamber and shallow rigid bodies deduced by high-resolution tomography studies.The authors support their observations by applying numerical modeling technique to simulate theoretical wavefront.Slowness anomalies observed at two arrays deployed on the same volcano were also described by Saccorotti et al. [2001].Similarly, Nisii et al. [2007] detected systematic differences between apparent slowness and backazimuth at Campi Flegrei volcanic area by using array data from an active source experiment.They attributed the results to a lateral velocity heterogeneity associated with the rim of the collapsed caldera structure.Slowness and backazimuth anomalies in teleseismic signals were also observed at Long Valley caldera [Steck and Prothero 1993].The largest backazimuth anomalies observed in volcanic areas are in the range between 50°, as found for Campi Flegrei [Nisii et al. 2007], up to 80° for Long Valley Caldera [Steck and Prothero 1993].
In this paper we analyze local and regional earthquakes recorded at Mt. Vesuvius by a seismic array temporary deployed in 1997-1998 and observe backazimuth anomalies up to 37°.We focus our attention on the backazimuth observed for the direct P-wave, estimated by the Zero Lag Cross-Correlation technique (ZLCC) [Frankel et al. 1991].The comparison of observed backazimuth with the theoretical direction of several local and regional earthquakes shows interesting differences that are investigated by computing the wavefront propagation through heterogeneous media.

Geological and geophysical setting
Mt. Vesuvius is a composite central volcano located in the Campania plain near the suburbs of Naples, Italy.It is formed by an ancient caldera (Mt.Somma) and a younger cone (Mt. Vesuvius).The eruptive history is characterized by long periods of quiescence alternated to periods of intense activity that usually start with strong explosive eruptions [Cioni et al. 1999, Santacroce et al. 2003, Cioni et al. 2008].Although it is quiescent since A.D. 1944, its activity is monitored in real time by modern networks and methods to mitigate the high volcanic hazard [Castellano et al. 2002, Baxter et al. 2008].In the last decades, geophysical investigation studies have been realized to improve the knowledge of the volcano structure.In particular, several tomographic studies aimed at determining the velocity structure have been realized [De Natale et al. 1998, De Matteis et al. 2000, Lomax et al. 2001, Scarpa et al. 2002].The common features evidenced by these studies consist of 1) a strong vertical and lateral P-wave velocity variations in the first 5 km below the crater [De Natale et al. 1998, Scarpa et al. 2002]; 2) the presence of high velocity core below the crater spreading towards southwest in depth [De Natale et al. 1998]; 3) the absence of a magma chamber in the first 4-5 km below the sea level.The complexity of the velocity structure and of the attenuation properties of Mt.Vesuvius area has been summarized by De Siena et al. [2009].The comparison between velocity and attenuation images shows the coincidence between low attenuation and high velocity zones in the central part of the volcano, at depth greater than about 0.5 km below sea level.The authors interpreted this result as solidified magma from the last eruption.

Data set
The seismic array was installed on the south-west flank of Mt.Vesuvius during a collaborative survey realized by Osservatorio Vesuviano and University of Granada, Spain [Del Pezzo et al. 1999].It was located 1.8 km from the summit crater, at an average elevation of 600 m a.s.l.It was composed by 18 vertical component and two three-component sensors deployed in an area of about 500 m in diameter (Figure 1).The natural corner frequency of 4.5 Hz of the sensors (Mark Product L15B) was electronically extended down to 1 Hz by a feedback circuitry.The homemade data acquisition system consisted of three laptop computers, each acquiring 8 channels.Data were acquired at 200 sps with a dynamic range of 16 bit, and time synchronization

Lat
Lon Depth (km) Source location and some other parameters of the 7 earthquakes discussed in this paper.Hypocenter and magnitude were taken from the catalog, while distance and theoretical backazimuth were computed from epicenter and array site coordinates.Earthquakes are ordered with increasing theoretical backazimuth { T .The observed backazimuth { O was estimated by array analysis by considering the first 3 s of signals (average values are indicated in bold inside the range of observed values).The last two columns contain the difference range between observed backazimuth and theoretical backazimuth, and between synthetic backazimuth and theoretical beckazimuth, respectively.was ensured by the GPS signal.A STA/LTA algorithm was applied to the incoming digital data in order to save only few minutes of signals for each trigger.The array was operative from November 1997 until June 1998.
The data set collected during the acquisition period consist of more than 450 earthquakes, tens of artificial signals and some rockfalls occurred inside the crater.Most of the earthquakes were of volcano-tectonic type, small magnitude (M D ≤ 2) and located inside the volcanic edifice.Some local earthquakes located outside the volcano and some regional earthquakes were also recorded (1.9 ≤ M L ≤ 5.5) and their detailed analysis constitutes the subject of this paper.The earthquakes of this group characterized by high signal-to-noise ratio (SNR) are listed in Table 1 and discussed throughout the paper.The localization of these earthquakes are taken from the INGV bulletin (ftp://ftp.ingv.it/pro/bollet).One example of array recording is shown in Figure 2.

Array analysis: method and results
Array data were analyzed with the Zero Lag Cross-Correlation (ZLCC) technique [Frankel et al. 1991] in time domain, to estimate slowness and backazimuth of coherent signals.The method is based on the search of the maximum of cross-correlation function (C max ) on a slowness grid.The values of slowness (Sx, Sy) corresponding to C max give the apparent slowness vector, from which slowness S and backazimuth are computed.ZLCC technique has been applied to many different frequency bands of one octave bandwidth (center frequencies of 1, 2, 3, 4.5, 6, 9 and 12 Hz) on short sliding window with 50% overlapping.The window length was set to include at least two periods corresponding to the central frequency of analysis.The high number of stations (at least 17 for all earthquakes resumed in Table 1) and the isotropic configuration of the array, associated with the extension of about 0.5 km (Figure 1), give to this array a very good resolution power for signals in the 2-8 Hz frequency range [Rost and Thomas 2009].At frequency lower than 2 Hz the resolution can be lower because the signal wavelength may be much greater than the array extension and because the correlation of seismic noise may be comparable with that of the earthquake signal.The results of array analysis are characterized by high correlation for many seconds at the earthquake onset in most of the analyzed frequency bands.The results for three earthquakes are shown in Figures 3, 4 and 5, while the best estimate of the P-wave backazimuth for the selected seven earthquakes are resumed in Table 1.The theoretical backazimuth { T , estimated from the array and epicenter coordinates, and observed backazimuth { O , are also reported in Table 1.For the observed backazimuth we report a range that include all values corresponding to the well correlated windows for the first three seconds of earthquake signal.We selected the phases with a correlation higher than a suitable frequency-dependent threshold, as shown by larger symbol in Figures 3, 4, 5.
We will describe in detail the results of ZLCC analysis for three earthquakes: 19980326 1626 (Figure 3), 19980404 0707 (Figure 4), 19980110 1923 (Figure 5).In Figures 3, 4 and 5 one seismogram is shown at the top to give an idea of the signal shape and SNR.The following plots show correlation, backazimuth and slowness with different symbols and colors for each analyzed frequency band.Results corresponding to windows characterized by correlation higher than an arbitrary chosen threshold (depending on the frequency of analysis) are shown by bigger symbols.A continuous line in the backazimuth plot indicates the theoretical backazimuth direction.The theoretical backazimuth of the earthquake 19980326 1626 (Figure 3) is { T = −27°, but the values of observed backazimuth are considerable different from the expected direction, with most of values in the range from −15° to 35°.For this earthquake the observed backazimuth is clearly frequency dependent, and the difference with the theoretical direction is greater at lower frequency.The observed backazimuth is computed as the average among all values corresponding to the windows whose correlation is higher than a threshold (big symbols in the figures) and considering the first few seconds of signal (from 12 to 15 in Figure 3) and all analyzed frequencies.The values of slowness are in the range between 0.2 s/km and 0.3 s/km.The results found for the earthquake 19980326 1626 are similar to the other earthquakes with theoretical WAVEFRONT DISTORTION AT MT.VESUVIUS backazimuth smaller than about 40 degrees, resumed in the first three rows of Table 1 (hereafter said group A earthquakes).
The earthquake 19980404 0707 (Figure 4) has { T = 51°.In this case only few windows at the signal onset are well correlated, and the corresponding values of backazimuth are very close to the theoretical value.The earthquake 19980110 1923 (Figure 5) has { T = 110°, but the ZLCC analysis shows observed backazimuth { O in the range from 80° to 105°, that is all values are smaller than expected.The values of slowness for well correlated windows are all very similar to each other, in the range between 0.2 and 0.3 s/km.The results found for earthquake 19980110 1923 are similar to the other earthquakes with theoretical backazimuth greater than about 70 degrees, resumed in the last three rows of Table 1 (hereafter said group B earthquakes).
The well correlated P-waves are characterized by slowness values between 0.2 s/km and 0.3 s/km for all analyzed earthquakes.The uncertainty on the slowness estimated by array methods depends mainly on the array extension, number of stations and SNR [Rost and Thomas 2002].By comparison of our array with other arrays scribed in literature [La Rocca et al. 2008], we believe that 0.03 s/km is a reasonable value for the slowness uncertainty.Considering the relationship between slowness uncertainty and backazimuth error [La Rocca et al. 2008], a slowness uncertainty of 0.03 s/km corresponds to errors of about 10 degrees on the observed backazimuths.However, looking at the results shown in Table 1, it is clear that for most of the analyzed earthquakes the difference between observed and expected backazimuth is much greater than 10 degrees.It is also noteworthy that the difference between observed and theoretical backazimuth is not systematic, being positive for group A earthquakes and negative for group B earthquakes (Table 1).Our experimental results for the 7 analyzed earthquakes can be synthesized in the following statements: 1) for the three earthquakes with theoretical backazimuth smaller than 40° we observe values greater than expected; 2) for the three earthquakes with theoretical backazimuth greater than 70° we observe values smaller than expected; 3) for the one earthquake with theoretical backazimuth of 51° the difference between observed and expected values is negligible.In this case the array site, the Vesuvius crater and the epicenter are aligned.
Out of the seven analyzed earthquakes, five are located on the opposite side of the volcano with respect to the array location (Figure 1; Table 1).For these earthquakes the raypath crosses a significant part of the volcano edifice before reaching the array.Therefore we Figure 4. Results of ZLCC analysis (correlation, backazimuth and slowness) for the earthquake 19980404 0707.For this earthquake the correlation is greater the the chosen thresholds only for about one second at the beginning.Note that the correlation appears shifted to the left with respect to the seismic signal due to the sliding window length.hypothesize that the main contribute to the observed backazimuth deflection must be produced by the propagation through the volcano structure.

Computation of synthetic wavefront in heterogeneous media
The investigation of the possible causes of the observed backazimuth deviations involves the simulation of wavefront propagation through a medium with well known characteristics, such as velocity discontinuities or velocity gradients.Wavefront distortion observed in other regions have been analyzed by taking into account complex velocity models [Shulte-Pelkum et al. 2003] and/or structural discontinuities in the crust [Arlitt et al. 1999].
The structure of Mt.Vesuvius volcano have been well studied in recent years, and many studies were focused on velocity tomography.Some works evaluated detailed velocity models using active and/or passive seismic data [De Natale et al. 1998, De Matteis et al. Lomax et al. 2001, Scarpa et al. 2002].Results exclude the presence of a magma chamber in the first 4-5 km below sea level, but significant heterogeneity in the velocity structure have been found by all authors.A common feature of all velocity models obtained by different studies is a positive velocity anomaly spanning the central part of the volcanic edifice.De Natale et al. [1998] identified a well defined trend in the NW-SE direction, with lower velocity in the northern sector and higher values in the southern sector.De Matteis et al. [2000] confirmed the presence of a high velocity body located NW of Mt. Somma-Vesuvius, and Scarpa et al. [2002] found strong vertical and lateral P-wave velocity variations in the first 5 km below the crater.In this last work two high velocity bodies can be roughly sketched: the first is located below the crater, and it extends from the sea level down to a depth of about 1 km [Scarpa et al. 2002]; the second one is located in the southern direction, possibly indicating a high concentration of dykes.All works are in agreement with the interpretation of a compact high density body located in the middle of volcanic complex [Berrino and Camacho 2008] and related to the oldest consolidated lava body.
By considering these results, we hypothesize that the positive velocity anomaly located beneath the crater can be the cause of wavefront distortion for the earthquakes observed at the array.Therefore we compute the seismic wavefront for different propagation azimuth in a medium with a positive velocity anomaly centered below the volcano crater.For simplicity we restrict our analysis to the horizontal plane, that is we assume a 2D velocity model and compute the propagation of the wavefront neglecting the vertical component.The result of such simulation, given the source epicenter, is a set of curves representing the wavefront at different travel times.Then the simulated backazimuth is computed as the direction perpendicular to the wavefront and directed toward the source.Synthetic wavefronts have been calculated by applying a finite difference approach to compute the horizontal propagation of seismic waves.Several methods have been developed to evaluate this topic [Vidale 1988, Ohminato andChouet 1997].We applied the finite difference method proposed by Vidale [1988] because it is easily implemented and requires low computing power.This method is developed for a velocity or slowness structure in two dimensions.The propagation of a wavefront, that is a curve with t(x,y) = constant, is computed from the eikonal equation that relates the travel time gradient to the slowness structure S(x,y).Computation is done on a regular 2D grid with equal step in the two directions.The grid step, h = 0.1 km in our analysis, is chosen from a trade off between computing time and resolution of the results.Our model is circularly symmetric, centered at the crater, and with slowness values starting from 0.2 s/km at the center and linearly increasing to 0.25 s/km at a radial distance of 4 km.At distance greater than 4 km from the center the slowness is constant, as shown in Figure 6.These slownesses are compatible with P wave velocity found by the different tomographic studies for the shallowest 5 km of crust beneath the volcano [De Natale et al. 1998, De Matteis et al. 2000, Scarpa et al. 2002].
Figures 7 and 8 show the results of wavefront simulation for each analyzed case.The position of the array is shown by a blue square, while the volcano crater is represented by a black circle.Black dashed lines represent the wavefront propagating through the homogeneous model, computed as circles centered at the source.Colored lines show the synthetic wavefronts computed by the finite-difference method with the non-homogeneous model.The wavefront distortion produced by the slowness anomaly is evident in all plots shown in Figures 7 and 8. Black arrows represent the theoretical backazimuth { T computed for the homogeneous model, therefore they always point exactly to the source.On the contrary, red arrows represent the synthetic backazimuth { S computed from the wavefront propagated through the heterogeneous model.In the region where the wavefront is modified by the slowness anomaly, the corresponding backazimuth is obviously deflected from the source direction, being perpendicular to the wavefront.The observed backazimuth { O , estimated by array analysis, is shown by the blue arrow, The difference be-tween { T and { S at the array site is evident in all analyzed cases, reaches values greater than 10° for some earthquakes.The lowest differences are observed for backazimuth of 51° and 35° (19980404 0707 and 19980422 2339 earthquakes in Figure 7).The difference <{ O > − { T between the mean observed and theoretical backazimuth, and { S − { T between synthetic and theoretical backazimuth are reported in the last two columns of Table 1.The sign of the difference { S − { T is in agreement with the results found from array analysis ({ O in Table 1).In Figure 8d the observed and synthetic backazimuth are plotted versus the theoretical backazimuth for all the seven earthquakes studied in this work.It is noteworthy as the observed values differ from the theoretical values more than the synthetic direction for all earthquakes of group A and group B. This is evident in the plots of Figures 7 and 8, where the red arrow (synthetic direction) is always located between blue (observed direction) and black (theoretical direction) arrows.The origin of X and Y axes is set at 40.821, 14.426, the center of Vesuvius crater, also shown by the black circle.The blue square indicates the position of the array.The four panels refers to the first four earthquakes of Table 1.

Discussion and conclusions
Array analysis has evidenced clear azimuthal anomalies for earthquakes spanning a significant range of backazimuth and epicentral distance.The backazimuth anomalies are in most of cases larger than both measurement errors (10°) and the variations associated with different time intervals and frequency bands.The possible causes of such discrepancies could be due either to a local heterogeneity beneath the array (local topography or structural discontinuities) and/or a velocity interface along the raypath.The observations that (a) azimuthal deflections are evident even at lowest frequency band, for which the correspondent wavelength are larger than array size and (b) the backazimuth deflections observed is not systematic but strongly dependent on the epicenter direction, suggest that the raypath bending occurs in a volume larger than the array area.The observed slowness values in the range 0.2-0.3s/km are normal for direct P waves of crustal earthquakes recorded at local and near regional distances.For farthest earthquakes, we expected to sample Pn phases characterized by horizontal slowness lower than 0.2 s/km.On the contrary, the results for the farthest earthquakes did not show differences in slowness with respect to closest events (0.2-0.3 s/km).This result can be another effect of the positive velocity anomaly which bends downward the upgoing incident rays (compared with the surrounding slower medium), with a consequent increase of the observed slowness.
Considering an average P-wave velocity of 3 km/s in the shallow layers below the array site, the range of incidence angle corresponding to 0.2-0.3s/km slowness is between 37° and 64°.This means that direct P waves of the earthquakes we have analyzed travel throughout the volcano structure in a wide range of depth, as shown in Figure 9. Assuming this range of incidence angle and considering some bending of the raypath produced by a velocity increase with depth, we have a rough estimate of the volume sampled by the studied earthquakes.Therefore we believe the model we have used to compute synthetic signals applies to the upper part of the volcano structure, roughly between 0 and 2 km depth below the crater.
The synthetic wavefront evaluated by finite differ- ence method confirmed how the backazimuth direction can be strongly distorted by a velocity anomaly along the raypath and represents a reasonable explanation of experimental data.It is clear that the sign of backazimuth deflection depends by three factors: 1) the array location with respect to the velocity anomaly; 2) the shape, size and other features (gradient velocity) of the velocity anomaly; and 3) the position of the seismic source.By considering the results of our analysis, we have shown the important role of velocity heterogeneity on wave propagation in the volcanic area of Mt.Vesuvius.More quantitative analysis should be carried out, such as the dependence of results from the frequency band (observed in several of the analyzed earthquakes), slowness and depth of seismic source.However, the small number of recorded earthquakes appropriate for the analysis described in the paper is the limit of our work.A quantitative analysis of the structural heterogeneity requires hundreds of earthquakes coming from a broad range of azimuth and distance.Considering the regional seismicity rate in southern Italy, the collection of such a database could require many years of data acquisition by an array like that installed in 1997.The availability of a rich data set would permit much more detailed analysis by taking into account the observed slowness and signal wavelength, and computing 3D simulation of the wavefront.Aware of the limit imposed by the small data set, our purpose was only a qualitative analysis of the observed backazimuth deflection.Deep investigation of the velocity anomaly, such as its shape, position, size, velocity contrast, and so on, is far beyond the scope of this work.Nonetheless, our qualitative results and a simple 2D simulation demonstrate once again the importance of seismic arrays for the study of local structures.Often small-aperture seismic arrays are used to locate the source of seismic signals by back-projecting the observed propagation direction of correlated signals.This is a typical use of arrays in the location of volcanic tremor or in the location of coherent noise sources.Our results suggest that much care must be given to such locations, taking into account that the signal wavefront may be considerable distorted by local heterogeneity.Therefore, the joint use of seismic arrays and a local seismic network, when available, is recommended for the location of seismic sources in volcanic environment.
In this work we have shown the important role of velocity heterogeneity on the wave propagation in the volcanic area of Mt.Vesuvius and the results are particularly relevant under the perspective of using array methods for monitoring purposes.

Figure 1 .
Figure 1.Map of southern Italy showing the position of Mt.Vesuvius (blue square) and earthquake epicenters (red star).The two insets show the topography of Mt.Vesuvius with the array location (a) and the array configuration (b).

Figure 2 .
Figure 2. Seismograms of the earthquake 19980422 2339 recorded at the array stations.The similarity of the waveform among the stations is evident only at the beginning and at lower frequency.Few seconds after the onset many seismograms look quite different each other, both in shape and amplitude.

Figure 3 .
Figure 3. Results of ZLCC analysis (correlation, backazimuth and slowness) for the earthquake 19980326 1626 obtained for six different frequency bands, as shown by different symbols and colors.At the top an unfiltered seismogram recorded by station 1D is shown.The results with correlation higher than a frequency dependent threshold are shown with larger symbols.The blue line in the backazimuth plot represent the theoretical backazimuth { T .

Figure 6 .
Figure 6.The slowness model used to compute synthetic wavefronts is characterized by a circularly symmetric anomaly centered at the volcano's crater (black circle).The blue square represents the array position.The origin of the reference system corresponds to the crater center (40.821, 14.426).

Figure 7 .
Figure 7. Synthetic wavefronts computed with the slowness anomaly (colored lines) and with a homogeneous slowness model (dashed black lines).Black and red arrows indicate backazimuth values obtained in the case of homogenous and non-homogenous slowness model.Bold bigger arrows show three different backazimuth at the array site: theoretical (black arrow), synthetic (red arrow) and observed (blue arrow).The origin of X and Y axes is set at 40.821, 14.426, the center of Vesuvius crater, also shown by the black circle.The blue square indicates the position of the array.The four panels refers to the first four earthquakes of Table1.

Figure 8 .
Figure 8. Plots a, b, and c show the synthetic wavefronts elaborated for the last three earthquakes listed in Table 1, similar to those shown in Figure 7. Plot d display observed and synthetic backazimuth versus the theoretical backazimuth for the seven earthquakes analyzed in this paper.

Figure 9 .
Figure 9. Cross-section of Mt.Vesuvius along the direction arraycrater.The two straight black line, corresponding to incidence angles of 37 and 64 degrees, give a rough boundary of the volume sampled by the analyzed earthquakes.Considering an increase of the wave velocity with depth, more realistic limits are suggested by the red dashed lines.