P-wave scattering and the distribution of heterogeneity around Etna volcano

Volcanoes and fault zones are areas of increased heterogeneity in the Earth crust that leads to strong scattering of seismic waves. For the understanding of the volcanic structure and the role of attenuation and scattering processes it is important to investigate the distribution of heterogeneity. We used the signals of air-gun shots to investigate the distribution of heterogeneity around Mount Etna. We devise a new methodology that is based on the coda energy ratio which we define as the ratio between the energy of the direct P-wave and the energy in a later coda window. This is based on the basic assumption that scattering caused by heterogeneity removes energy from the direct P-waves. We show that measurements of the energy ratio are stable with respect to changes of the details of the time windows definitions. As an independent proxy of the scattering strength along the ray path we measure the peak delay time of the direct P-wave. The peak delay time is well correlated with the coda energy ratio. We project the observation in the directions of the incident rays at the stations. Most notably is an area with increased wave scattering in the volcano and east of it. The strong heterogeneity found supports earlier observations and confirms the possibility to use P-wave sources for the determination of scattering properties. We interpret the extension of the highly heterogeneous zone towards the east as a potential signature of inelastic deformation processes induced by the eastward sliding of flank of the volcano.


Introduction
The TOMO-ETNA field project is a joint research effort of theInstituto Andaluz de Geofísica -University of Granada, Spain, the Istituto Nazionale di Geofisica e Vulcanologia (INGV), Italy, and the GFZ German Research Centre for Geosciences in Potsdam, Germany [see Ibáñez et al. 2016aIbáñez et al. , 2016b, in this volume], in this volume].The project is part of the European Research Project MED-SUV (Mediterranean Supersite Volcanoes).Aim of the terrestrial and marine experiment TOMO-ETNA is to resolve the structure beneath the Etna volcano and the surrounding areas of Sicily in a never reached high structural resolution using active and passive seismic methods.In the frame of the onshore measurements seismometer stations were deployed around Etna and across Sicily.The recording instruments were provided by the Geophysical Instrument Pool Potsdam (GIPP) established at the GFZ Potsdam in Germany [Ibáñez et al. 2016a[Ibáñez et al. , 2016b, in this volume], in this volume].
We will fist introduce the methodology, the geological and seismological setting of this study at the Etna volcano, and observations at other volcanoes in Sections 1.1 -1.3.In Section 2 we present our data selection.In Section 4 we will show our results of the different scattering properties around Mount Etna.The results will be discussed in Section 5. Conclusions are given in Section 6.

Methodology
Besides the large scale seismic velocity and attenuation structure the small scale seismic velocity structure is also of great interest especially in volcanic settings.Small scale here refers to a structural length scale of the order of typical seismic wavelength or below.Details at such resolution cannot be inferred deterministically with seismic methods but related structural in homogeneities strongly influence the propagation of seismic waves by scattering, i.e. a redistribution of seismic energy.This redistribution of energy is the more efficient the more small scale heterogeneity is present in the subsurface, which is responsible for the scattering processes.In volcanic environments the presence of magmatic intrusions, lava flows and ash layers and other complexities in the internal structure usually cause strong wave scattering in volcanic settings.
Wave scattering has different effects on seismic waveforms that are discussed in a comprehensive review by Sato et al. [2012].The most prominent effect is the generation of the seismic coda that was first discussed by Aki and Chouet [1975].They demonstrated that the seismic energy that arrives at a seismic station significantly after the arrival of ballistic P-and S-waves thereby forming the seismic coda is best explained by waves that are scattered at distributed inhomogeneities in the medium.This energy is removed from the ballistic seismic wave field just like ballistic energy is decreased by intrinsic attenuation.Distinguishing between intrinsic and scattering attenuation therefore is not possible if only the amplitude decay of direct waves is considered [Sato et al. 2012].This complexity illustrates the importance of coda-wave analysis.Another consequence of wave scattering is envelope broadening.This effect refers to the progressive broadening of the direct P-and S-wave pulses due to perturbations of the ray direction caused by forward scattering.Sato [1989] studied the broadening of seismogram envelopes in relation to the random inhomogeneity distribution in the lithosphere beneath Japan.Saito et al. [2005] explored the envelope broadening and the maximum-amplitude decay of high-frequency seismograms in a volcanic area in Japan by modeling the observations with the Markov approximation that accounts for multiple scattering in the forward direction.An important characteristic of envelope broadening is the peak delay time that describes the time difference between the first arrival of a certain phase and the arrival time of its energy maximum.Obara and Sato [1995] initially observed that the S-wave peak delay strongly increases in Japan when the waves cross the volcanic front in northern Honshu.The precise location of the strong scattering bodies was found later by Takahashi et al. [2007] who could show that only waves that propagate directly through the Quartenary volcanoes experience strong envelope broadening.Calvet et al. [2013] measured peak delay times of around 10,000 S-arrivals of earthquakes in the Pyrenees.The large number of observations allowed to map the distance corrected peak delay times as a proxy of the spatial distribution of scattering strength.This revealed significant spatial differences in the small scale heterogeneity that correlate with geological units but were not recovered in the large scale seismic velocity structure.
In our study we investigate spatial differences in the small scale heterogeneity distribution around Etna volcano using the records of air-gun shots in the Ionian and Tyrrhenian seas performed during the TOMO-ETNA project.In contrast to earthquakes that emit seismic energy predominantly as S-waves with a clear radiation pattern, the air-guns are pure P-wave sources with isotropic radiation.Because of these source characteristics we analyze the P-wave phase.The different influences of wave scattering on P-and S-waves require a new approach to the analysis of the waveforms.In analogy to the mapping of peak delay times used by Calvet et al. [2013], we follow a purely qualitative approach that is based on the details of the scattering process.An important property of elastic scattering is that conversion from compressional waves into shear waves is much more efficient than scattering from shear waves into compressional waves.This effect influences the equipartition ratio of the scattered wavefield -a fundamental property that describes the ratio of the energy in the two wave modes with a strong dominance of the shear mode.
Whereas the S-wave envelope broadening observed as peak delay is caused by multiple scattering of the large amount of S-waves in the forward direction, the P-wave pulse is mainly attenuated by scattering as energy -if scattered -is predominantly converted to S-waves.Due to the slower S-wave velocity, the scattered phases arrive later in the recorded waveforms.If these scattered waves propagate along the shortest path between source and receiver, the scattered energy arrives within a time window between the P-arrival, if the scattering occurred close to the receiver, or shortly before the S-arrival, if scattering is close to the source.This relocation of seismic energy from the direct Ppulse into these later parts of the seismogram is a sensitive measure of the scattering strength.We quantify this process as the ratio between the P-pulse amplitude and the averaged amplitude of the seismogram envelope around the arrival of the direct S-waves.This approach has some similarity with the multiple lapse time window analysis (MLTWA) proposed by Fehler et al. [1992].This technique uses the ratios of the seismic energy in three S-wave coda time windows and allows to for quantitative estimates of the intrinsic and scattering attenuation based on extensive numerical model calculations [e.g.Giampiccolo et al. 2006;Carcolé and Sato 2010].Due to the strong dissipation and weak impulsive earthquake sources inside volcanoes MLTWA cannot be applied at volcanoes in its original form.A modeling of the S-coda together with the amplitude of the direct S-wave was used by Sens-Schönfelder and Wegler [2006] to estimate scattering and attenuation parameters in continental crust.
The use of high frequency P-wave coda signals has mostly been restricted to teleseismic signals [e.g.Korn 1988;Rothert and Ritter 2000;Gaebler et al. 2015].In this study we propose to use P-waves and their coda at source distances smaller than 100 km.In this first attempt we explore the regional variations of the scattering effects rather than quantifying them with an absolute reference.
The following questions will be addressed: (1) how far can methods developed for the analysis of scattered waves from earthquakes be transferred to the air-gun shots and (2) are there regional differences in the scattering properties at Etna volcano.
To characterize these variations of the attenuation mechanisms underneath Etna volcano, we will use a dense seismic network, which consists of nearly 100 (?) temporary stations provided by the GFZ.

Seismic wave scattering at volcanoes
The internal structure of volcanoes results in strong scattering of seismic waves.This has been recognized at Etna and Campi Flegrei (Italy) by Del Pezzo et al. [1996] using the approximation of multiple S-wave scattering in a half space.The diffusion model that assumes strong scattering in a half space has been applied to data from Merapi volcano (Indonesia) by Wegler and Lühr [2001].Extensions of this model to more complex boundary conditions were provided by Wegler [2003] who used a strong heterogeneous layer above a weakly scattering half space.The same approach applied Prudencio et al. [2013aPrudencio et al. [ , 2013b] ] to different volcanic environments and by using active sources.Friedrich and Wegler [2005] modeled the volcano as a strongly heterogeneous cylinder embedded in a homogeneous half space.Topography has been incorporated in the modeling of wave scattering in Merapi volcano by Parsiegla and Wegler [2008].Without providing detailed information about the distribution of heterogeneity in the vicinity of volcanoes the investigations clearly showed that scattering of seismic waves in the volcanic edifice is far strongly than in the surrounding crust.For a comprehensive review of these investigations we refer to the review paper Del Pezzo [2008].
More detailed information about the distribution of heterogeneity inside and around volcanoes result from tomographic approaches as presented by De Siena et al. [2010] for Campi Flegrei, Vesuvius [De Siena et al. 2014a] andMount St. Helens [De Siena et al. 2014b].With a similar approach and by using active sources Prudencio et al. [2015aPrudencio et al. [ , 2015b] ] imaged the distribution of heterogeneity of Deception Island and Tenerife (Canary Islands).A comparable analysis is not available yet for Mount Etna.

Structure of the Etna volcano
Mount Etna is an active stratovolcano with about 3329 m height, which makes it the highest active vol-cano in Europe.The volcano is located at the east coast of Sicily (Italy), close to the city of Catania which is endangered by its volcanic hazards.Because of the high frequent activity phases, Mt Etna is a useful object to study volcanic processes.In the following, we describe the main geological and seismological structures of the area.
Mount Etna is located in the middle of an intersection of two active master faults [Azzaro 1999;Acocella et al. 2013].First there is the "Malta Escarpment" (Figure 1), which is the boundary between the continental lithosphere of Sicily and the oceanic lithosphere of the Ionian Sea.Second there is the "Messina-Fiumefreddo" line, which is striking NNW and NE, and which is mainly responsible for the location of the Etna volcano.Also, there is the "Apennine -Maghrebian" thrust belt.Because of the interaction of the regional tectonics and its stress field, the area around the volcano is influenced by many local earthquakes.
Especially the tectonic settings of Etna's flanks are important for observations.In the eastern part of the volcano there is a clustering of tectonic lines.These lines reflect a north oriented prolongation of the "Malta Escarpement", and imply an instability of the eastern flank of the volcano.Because of the interaction between these flank instability and gravity forces there is a seaward sliding of the volcano's eastern side [Azzaro 1999].The near surface dynamics of Mount Etna is mainly controlled by the flank instability.

Data acquisition
For the TOMO-ETNA campaign the GIPP provided 77 3-D geophones with 4.5 Hz natural frequency plus 3 Mark L4-3D 1 Hz short period seismometers in combination with DATA-CUBE3 data-loggers, and 20 Trillium Compact broadband seismometers (ground velocity sensors from 120 sec to 100 Hz) in combination with EarthData PR6-24 recorders.The sample rate of the instruments was set to 100 samples per second.The network was supplemented by 70 permanent stations operated by INGV, leading to a total number of 168 seismic stations onshore.In addition, 15 ocean bottom stations (OBS) of the Spanish UTM-CSIC research centre as well as 12 Italian OBS were deployed offshore at the bottom of the Tyrrhenian Sea und Ionian Sea by the research vessels "Sarmiento de Gamboa" and "Galatea".95% of all onshore instruments were recording in the time between June 18 and June 26, 2014.During the experimental phase 11 of the 80 short period stations were relocated and in total data from 91 recordings sites could be used for this study (Figure 1).
The active seismic experiment started offshore on June 27 at 17:25 h GMT with air-gun shooting from the research vessel "Sarmiento de Gamboa" with a series of low intensity shoots in intervals of 30 seconds.The air-gun shots ended on July 16.After this date the 80 short period stations were removed from the field.In total 13,000 air-gun shots were recorded.After the active experiment phase only the deployed broadband stations remained operative in the field until end of October 2014.

Data pre-processing
The data were available in an instrument specific format.Because of the usage of normalized envelopes in a small frequency range, a deconvolution of the station response was not necessary.An instrument correction is therefore obsolete for this specific use.
First, we extracted the air-gun shots from the data according to the shooting times.We tapered the ends of the shot windows with a cosine taper and used a fourth-order Butterworth bandpass filter from 5 Hz to 10 Hz to select the most energetic frequency band.
Afterwards we calculated the envelopes by means of the Hilbert transformation as the amplitude of the analytic signal.Finally we subtracted the noise level defined as the average envelope amplitude of the first 5-10 percent of the recording trace prior to the P-wave arrival from the envelopes.

Data quality
In contrast to stations located at the east coast of Sicily close to the air-gun shots, stations located around Mount Etna's summit have a poor signal-to-noise ratio on all components.Clear onsets of different shots were not identifiable and furthermore, several recording traces contain signals of unclear origin.Probably in addition to the signals of the air-gun shots, the traces contain tremor signals or local earthquakes as well as anthropogenic noise.A large part of this work was therefore devoted to separating shot signals from disturbing signals in order to improve the signal-to-noise ratio at the stations with a large distance to the shot positions.
To improve the signal to noise ratio after the preprocessing, we wanted to use a large number of shots for stacking.We therefore defined areas to group closely located shots for stacking.These areas comprise ring segments defined by a maximum distance variation of ±5% and an azimuthal variation of maximum 4 degrees between the recorder site and the shot positions.Thus the ring segments grow with increasing distance which allows us to stack a larger number of shots for larger distances.The ring segments are located around the intersection points of the shot lines to maximize the shot density in the stacking areas as illustrated in Figure 1.
A first visual inspection of the wave forms showed that the air-gun shots did not generate enough S-wave energy to actually study S-wave codas.This is a conse- quence of the small amount of energy released by the air-gun sources that is restricted to compressional waves.Therefore an estimation of the quality factor of the coda decay Q c was not possible with our dataset.As can be seen in the envelope section of station E078 in Figure 2 there is no clear S-wave onset visible in the data.Even the P-wave arrival at this station behind the volcano is characterized by an emergent increase of energy rather than a sharp impulse.

Data processing
In the present study the principal observation is the amplitude ratio between the direct P-pulse and the following coda.We call this value "coda energy ratio".To estimate this ratio we could use the maximum of the P-wave train as proxy for the direct wave amplitude.However, due to wave interference and noise there is a significant scatter in the maximum amplitudes of the P-pulses that can be reduced by fitting a smooth curve to a short time window around the P-wave arrival.As we also lack precise arrival time information for the Pwaves, we developed a measurement procedure that simultaneously estimates the arrival time, the peak time and the peak amplitude.
To find the onset of the direct wave and to pick the energy maximum, the following empirical function was fit to the recordings: where t is the lapse time, x is the distance between source and the station, v is the estimated velocity of the P-wave, d is the peak delay time and H is the Heaviside function.A is an amplitude scaling factor.Function ( 1) is zero for t < x/v and the maximum of f = A is assumed at t = x/v + d.For large t the function approaches zero again.This function is fit in a least squares sense to the stacked envelopes in a grid search over a number of reasonable values for the velocity v, amplitude A, and the peak delay d.The setting with the smallest misfit was adopted as empirical fit to the recordings (Figure 3) resulting in a measurement of the direct P wave onset, its amplitude and the peak delay time.
In addition we defined time windows that divide the seismogram into a direct P-wave time window and a P-coda time window.Inside these time windows we estimate the mean value of the wave energy and calculate the ratio between the energy of the direct wave and the energy of the coda.This ratio is a proxy for seismic scattering.A high coda energy ratio indicates strong scattering along the seismic ray path.For a low ratio, a large amount of energy reaches the station directly after the P-wave onset, indicating only weak scattering.
For the definition of the time window length, we used an apparent wave velocity to include the distance dependence and to avoid overlapping of the different time windows.For the direct P-wave, we used a time window in which the start is the picked onset time of the fit-function (1) and the end is determined by an apparent propagation velocity of 3.3 km/s.The second window, which measures the energy of the coda, is determined by velocities of 2.0-3.0 km/s.The definition of the time window has a small impact on the measurements of the coda energy ratio.This point will be discussed in Section 5.
With increasing distance the values for the coda energy ratio are characterized by a decreasing trend.
, , , This is illustrated in Figure 4.The origin of this effect is the attenuation in the increasing time difference between the two time windows.We correct our observations for this background effect with a linear regression of the logarithmic ratio over distance (Figure 4).

Results
First, we illustrate the results of the coda energy ratio (P-coda/P-onset) by plotting them over the respective station site (Figure 5).The stations are ordered from east to west across the volcano.The different values for one station site indicate the dependence of the coda energy ratio on the back azimuth of the incident rays in 45 degree intervals.Most of the values are widely spread, but there is a tendency that the eastern most stations show generally higher values.This implies an area in the eastern part of the volcano with strong scattering.
To better illustrate the results of the energy ratio calculation, we assume that the sensitivity of the meas-urements is concentrated along the ray path.This is motivated by the following argument: The amplitude of the coda time window reflects waves that traveled for a long time in the medium capturing influences from a broad region.Especially late coda waves are not sensitive to the specific geometry of the source receiver configuration which is the basis of the coda normalization method [Aki 1980].In contrast the P-pulse amplitude can only be influenced by material properties along the direct ray path.The resulting sensitivity of the ratio thus reflects the ballistic propagation and is concentrated along the ray path.For this reason we plot the results at the position of the station sites depending on the back azimuth of the incident ray (Figure 6).We split the azimuth into eight bins of 45° aperture and indicate the measurement values as colored wedges pointing in the respective back azimuth directions.Results are only included, if at least five measurements are available for averaging in a 45° bin.Yellow-and redcolored pie wedges define a high ratio, green-respec-   To increase the confidence in these observations we compare the energy ratio with observations of the peak delay time of the P-phase.Similar to the energy ratio the peak delay is a qualitative measure of the scattering along the ray path.Figure 7 shows a scatterplot of the energy ratio versus peak delay.The correlation between the energy ratio and the peak delay time can be recognized by the linear trend of the least squares fit in Figure 7.

Discussion
The results given in the previous section are in accordance with other observations of scattering in stratovolcanoes [e.g.Wegler and Lühr 2001;Del Pezzo 2008;De Siena et al. 2010;Prudencio et al. 2013a].The heterogeneity of the volcanic edifice is reflected by the high coda energy ratio observed on and in terms of the seismic ray paths behind the volcano.This heterogeneity is most probably caused by the alternated deposition of intrusions, lava and ash layers.
In the northern part of the volcano the energy ratio is very low (green-colored in Figure 6).Here the direct P-waves carry the main energy and sparsely scattered energy reached the stations, which indicates the absence of strong heterogeneity in this area.
In contrast the stations located east of the volcano that have almost ideal noise conditions in the close proximity to the sources, are characterized by high coda energy ratio, indicating strong subsurface heterogeneity at the east coast of Sicily.In Figure 8 we mark the area that is affected by elevated scattering.There is no precise information about the depth range in which this heterogeneity is located.Based on the simple argumentation that the sensitivity follows the direct ray path plus the concentration of sensitivity at the station we can expect that the observations are mainly related to the structure of the upper crust dominated by the shallow part.
We think that the increased heterogeneity that causes this observation might be related to the seaward sliding of the volcanic rock mass [Azzaro 1999].
The methodology used here is based on the assumption of a homogeneous background velocity in the subsurface.Deterministic structure, e.g.large-scale seismic velocity anomalies, can thus result in a systematic bias of the results.A potential source of such a bias are reflected waves such as the PmP-phase which is reflected off the Moho discontinuity.If this phase adds significant seismic energy in one of the time windows,  used to calculate the energy ratio, it could influence our observations.The depth of the Moho underneath Sicily is about 25-30 km [Grad et al. 2009].Like Gajewski and Prodehl [1987] showed for similar Moho depths in the Black Forest, the PmP-phase arrives at a distance of 40 km about 4 s after the direct P-wave.Accordingly this phase can theoretically influence our coda window.We used two approaches to estimate the PmP influence on our measurements.At first we visually checked the record sections for phases propagating with PmP-velocity.We could not identify a notable energy increase at the theoretical PmP arrival times.As a second approach to test PmP interference, we restricted observations of the energy ratios to a distance of less than 40 km (Figure A1).At these distances, the PmP-phase should have no effect on our results, because it arrives after the time windows for the analysis.In fact, there is no clear difference between these results and the results obtained from all available distances.For this reason we can exclude an impact of the PmP-phase on our results.
To test the impact of the time window choice, we performed a parametric study.In particular we used shorter time windows for the direct wave and larger time windows for the coda.For the following time windows we obtained results shown in Figure A2.The end of the new time window for the direct P-wave is now determined by an apparent propagation velocity of 4.0 km/s.The second time window for the energy of the coda is now determined by velocities of 3.8-2.0km/s.The spatial distribution of the ratios is nearly the same as the results presented in Section 4 (Figure 6).We also tested other choices of time windows and we conclude that the method is stable with respect the time window selection provided that the windows are chosen such that they contain signal energy above the noise level.

Conclusion
Motivated by the observation of significant lateral variations of scattering properties in the Pyrenees by Calvet et al. [2013], we analyzed recordings of air-gun shots around Sicily to infer the distribution subsurface heterogeneity in the vicinity of Etna volcano.As the sources are purely compressive pulses, we cannot analyze the peak delay of S-waves or measure the quality factor of the late time S-wave coda decay.Instead we used stacks of air-gun shot records to measure the energy ratio between the direct P-phase and a later time window in the P-coda.These observations are interpreted in terms of the scattering strength along the ray path of the direct wave.
We find that the coda energy ratio is stable with respect to changes of the choice of the coda time window.It correlates with the peak delay of the P-phase which is an independent measure of the scattering strength.
In our results Mount Etna is marked as a pronounced region of increased heterogeneity that strongly scatters seismic waves.Towards the north of Etna volcano we find a more homogeneous crust with relatively weak scattering.This is in contrast with the regions east and partially south-east of Mount Etna that are characterized by high coda energy ratios indicating increased heterogeneity.This area is characterized by local complicated geological structures and our results indicate the presence of strong inhomogeneities, potentially related to the seaward sliding of the eastern flank of the volcano.
Our investigation demonstrates that the signals of air-gun shots can be used to investigate subsurface scattering properties.Even though we cannot use measures like S-wave peak delay or S-wave coda Q due to the low generation of S-waves, we can make use of the P-phases to infer the transfer of energy form the ballistic P-phase towards later parts of the seismogram.This transfer can only occur due to scattering.A proper mapping of the scattering strength with a quantitative tomographic inversion is not possible yet, because there is no clear formulation for the sensitivity kernels of the elastic scattering process available.This is therefore subject of future work.

Figure 1 .
Figure 1.Overview of the main geological structure of Mt.Etna, as well as the used station network (red triangle) and the intersection points used for locating the ring segments for stacking (yellow stars) on the shot lines (black lines).Data examples are shown for station E078 (blue triangle).

Figure 2 .
Figure 2. Sumtraces of envelopes over distance.Shown is the horizontal component (N-S) of station E078.The traces are normalized relative to its respective maximum.The red line symbolize the onset of a P-wave with estimated velocity of v P = 4.4 km/s.There is no clear S-wave onset or coda.

Figure 3 .
Figure 3.At the top is a data example of a stacked trace (blue) and the fitted function (red) from station E078.The two horizontal lines are the meanvalues of the different time windows.At the bottom you can see the single traces which were used for stacking.

Figure 4 .
Figure 4. Logarithmic energy ratio plotted over distance.The red points symbolize the non-corrected mean values with a decreasing trend, which is most probably a background effect.This effect was corrected with a linear regression (blue stars).

Figure 5 .
Figure 5.The energy ratio plotted over their respective station from east to west of the volcano.Different values on one station symbolize the mean value of a 45 degree interval.

Figure 6 .Figure 7 .
Figure 6.Energy ratio, depending on back azimuth.Each pie wedge uses at least five single values for averaging.Green and blue colors symbolize a low energy ratio (weak scattering), yellow and red symbolize high energy ratio values.

Figure 8 .
Figure 8. Area of elevated scattering caused by increased heterogeneity.