Quantifying strong seismic propagation effects in the upper volcanic edifice using sensitivity kernels

In volcanic environments, the correct interpretation of the signals recorded by a seismic 8 station is critical for a determination of the internal state of the volcano. Those signals 9 collect information about both the seismic source and the properties of the path travelled 10 by the seismic wave. Therefore, understanding the path eﬀect is necessary for both source 11 inversions and geophysical investigation of the volcanoes’ properties at depth. We present an 12 application of the seismic adjoint methodology and sensitivity kernel analysis to investigate 13 seismic wave propagation eﬀects in the upper volcanic ediﬁce. We do this by performing sys-14 tematic numerical simulations to calculate synthetic seismograms in two-dimensional models 15 of Mount Etna, Italy, considering diﬀerent wave velocity properties. We investigate the re-16 lationship between diﬀerent portions of a seismogram and diﬀerent parts of the structural 17 volcano model. In particular, we examine the inﬂuence of known near-surface low-velocity 18 volcanic structure on the recorded seismic signals. Results improve our ability to understand 19 path eﬀects highlighting the importance of the shallowest velocity structure in shaping the 20 recorded seismograms and support recent studies that show that, although long-period seis-21 mic events are commonly associated with magma movements in resonant conduits, these 22 events can be reproduced without the presence of ﬂuids. We conclude that ediﬁce hetero-23 geneities imparts key signatures on volcano seismic traces that must be considered when 24 investigating volcano seismic sources. 25


Introduction
Complex stratigraphy, plumbing systems filled with magmatic fluids, pronounced topography and a broad set of possible source scenarios can all be found at most volcanoes.In order to obtain information about the internal state of a volcano, volcano seismology tries to understand how the above elements are seen in the signals recorded at the seismic network.
In previous papers, Bean et al. (2008) accounted for near-surface velocity structure through the computation of synthetic seismograms and sensitivity kernels to help with the interpretation of source inversion results for LP seismic events at Mount Etna, Italy (Lokmer et al. (2007)).Here, we extend that work to investigate the relationship between different parts of a seismogram and different portions of the structural model.In order to achieve this, we use the software SPECFEM2D (Tromp et al. (2008)) to perform 2D numerical simulations of seismic wave propagation in the presence of topography and heterogeneous velocity structure.Then, we take advantage of the adjoint methodology presented in Tromp et al. (2005) to study the sensitivity of seismograms with respect to P and S wave velocities by calculating the associated travel time sensitivity kernels.These kernels will show the regions of the velocity model that mostly affect the wavefield arriving at the stations within the seismogram time window under investigation.This yields important insights into how path effects control the overall structure of the recorded seismograms.

Method
Using inversion and adjoint methodologies (e.g.Tarantola (1984), Tarantola (1987), Tarantola (1988), Talagrand and Courtier (1987), Crase et al. (1990), Akcelik et al. (2002), Akcelik et al. (2003)) and connecting them with finite-frequency traveltime tomography (e.g.Zhao et al. (2000), Dahlen and Baig (2002) and Hung et al. (2000)) Tromp et al. (2005) demonstrated that Fréchet derivatives of objective functions, minimizing differences between traveltimes, may be written by means of traveltime sensitivity kernels.These kernels are built from the wavefield for a given model and a wavefield obtained by using time reversed signals as simultaneous sources at the receivers.The methodology allows for a better understanding of the origin of all phases in a given seismogram with respect to the velocity structure through which waves propagate.This is important in the context of highly heterogeneous volcanoes, where long duration seismic arrivals are common.We start with a brief overview of methodology.
In an isotropic medium m defined by the density (ρ), the shear (µ) and bulk (κ) moduli, the displacement wavefield s(x, t) induced by a seismic source can be described by the equation of motion where x denotes the spatial coordinates, t the time, f the force representing the seismic source and σ the stress tensor satisfying the constitutive relationship with the strain tensor In these equations, ∂ 2 t indicates the second partial derivative with respect to time t, ∇ the divergence operator and I the identity matrix.
In order to minimize differences between the calculated data at N seismic stations x r , r = 1, . . ., N and the observed signals, Tromp et al. (2005) introduced the traveltime misfit function where T r (m) denotes the synthetic traveltimes recorded at x r for the model m and T obs r the observed traveltimes.The traveltime Fréchet derivatives of this misfit function, representing the sensitivity of a seismogram with respect to the model parameters, can be expressed as where V denotes the model volume, δ the relative model perturbations and Kρ , Kµ and Kκ are the so called banana-doughnut kernels given by with [0, T ] indicating the time interval under investigation, s † the traveltime adjoint field generated by time-reversing the predicted components of ground velocity at receiver r, D the strain deviator tensor and D † its traveltime adjoint.
Alternatively, a representation of the kernels in terms of density, compressional wave speed (α) and shear wave speed (β) may be written as determined by These derivatives may be obtained by means of two numerical simulations for each seismic source: one calculation for the current model (forward) and a second (adjoint) calculation that injects the time-reversed signals in the receiver under investigation.Construction of travel time sensitivity kernels for the P-SV wave arrivals is based upon the interaction between forward and adjoint wavefield defined by equations 5.
This methodology is accommodated in SPECFEM2D (Tromp et al. (2008)), a twodimensional elastic wave propagation code based upon the spectral-element method (SEM).
We use that implementation in the context of complex wave effect in heterogeneous volcanic environments (i) to determine the region of the structural model that dominates the contribution to different portions of the seismograms and (ii) we check the effects of different source depths source mechanisms and frequency content.

Models and data
We use the methodology outlined in the previous section to study wave propagation in the upper edifice model of Mount Etna, Italy (Figure 1a).The 2D model is 20 Km x 6.25 Km. Digital topography was obtained from the CGIAR consortium for spatial information  Christensen (1996), Cayol and Cornet (1998), Bonaccorso andDavis (1999), Gercek (2007), Salah and Seno (2008)), we fix it at 0.25 to obtain the corresponding P-wave velocity structure.While the Vp/Vs will affect the more subtle details of the seismograms, it will not influence the core finding outline here.
Finally, from the P-wave velocity model, we derived the density model using the Gardner's relationship (Gardner et al. (1974)) where v p is measured in m s −1 and density in kg m −3 .The heterogeneous velocity model is summarized in Figure 1b and shown in Figure 1a.In this paper we focus on low frequency seismicity as it is often seen to proceed volcanic eruptions.As source mechanisms we use a Ricker wavelet applied either as a single vertical force or as an explosion and we consider dominant frequencies of 0.7 or 1 Hz (Figure 1c).
Even if the existence of the single force in LP sources is debatable (De Barros et al. (2013), Saccorotti and Lokmer (2021)), it is reported as the part of MT inversions (e.g.Kumagai et al. (2002), Kumagai et al. (2005)), and we use it here as case scenario.On the other hand, while the central frequencies of 0.7 and 1 Hz sources are similar, we chose them to test the sensitivity to small central frequency differences in the most common range of LP frequencies.We focus on two depths.One where the source is shallow and in the low velocity layer (113 meters depth), the other where the source is deep and lies outside the low velocity layer (3613 meters depth).Locations are denoted by black stars in Figure 1a.

Results
We perform full wavefield numerical simulations calculating synthetic seismograms for each model, for two different source frequencies and depths.We also calculate the corresponding travel time sensitivity kernels for successive time windows sliding along our seismograms.
These kernels show the region of the velocity model that affects the wavefield which is arriving at a given station within a defined time window.
Here, we present the most significant results from our simulations.
4.1 LP case (0.7 Hz) First, we show results for simulations considering the source Ricker wavelet applied as a vertical force with dominant frequency of 0.7 Hz.Later we will compare the results of a vertical source with those of an explosion.Volcanoes often have a velocity structure comprising very slow near surface deposits and a strong velocity gradient in the top 500 m to 1 km (e.g.Ferrazzini et al. (1991), De Luca et al. (1997), Chouet et al. (1998), Mora et al. (2006), Cauchie and Saccorotti (2013)).Due to the fundamentals of wave propagation, such complex volcanic structures produce complex waveforms, which may lead to misinterpretations of processes that generate them.For instance, source resonance-like signals may instead be structure-induced reverberations.This can be also a problem for source inversions, where the propagation effect needs to be modelled by numerical simulations.This requires a detailed knowledge of the velocity structure, which is usually not available in sufficient detail to allow path effects to be properly removed, leading to contamination of the source solution.Here we are not making any inferences about the actual source mechanism, merely assessing the sensitivity to the source (see Bean et al. (2008) and Trovato et al. (2016) for discussion on sources).
For the source at shallow location (113 meters depth) applied as a vertical force, Figure 2 denotes the velocity seismograms corresponding to simulations for both a homogeneous (Figure 2a) and the heterogeneous (Figure 2b) model.We see, in the homogeneous case, that original signal signature is affected by topography only, while, in the heterogeneous case, we observe the strong effect of the velocity structure on the traces producing, in addition, an extended wave train in the seismograms.Note that the length of the seismograms increases with the distance from the source (Figure 2b).The same effect, on real data from Mt Etna, was observed by Bean et al. (2014).
In particular, for the signal recorded at position -6000 meters (first seismogram from the left), we plot in Figure 3 the corresponding velocities and associated spectral content and spectrograms showing that, for the layered case, the signal has a higher frequency onset followed by a harmonic waveform containing one or up to several dominant periods.Also, Figure 3 shows a considerable increase in signal amplitude when near surface layering is considered, an increase which is more pronounced for the horizontal component.For the same station, in Figure 4 we display the corresponding P-and S-wave sensitivity kernels for consecutive 3-seconds time windows of the seismograms, showing that, for the homogeneous case (Figure 4a), the first seconds are mainly controlled by the area between the source and the receiver yielding the typical banana-doughnut shapes (Tromp et al. (2005)) distorted by the topography effects.In the layered case (Figure 4b), the sensitivity area remains within the low-velocity layers.This is a key observation as it demonstrates that short duration shallow sources can produce long seismograms that are almost entirely a consequence of the details of near surface velocity structure.
For the signal recorded on the highest station (eighth receiver from the left) Figure 5 shows the corresponding velocities and associated spectral content and spectrograms.
We observe that for the vertical component, even for the layered case, the source signature prevails with only small path effects.However, as seen also in Figure 3, signals are amplified due to the low velocity layers.Specifically, we observe the notable increase of the role played by the horizontal component in the layered case.We observe that, for the homogeneous model the mean vertical component increases by a factor of 2 if a near surface layering is considered, while the horizontal component increases by a factor of 10.

Sensitivity to source depth
For the case of the 0.7 Hz Ricker wavelet applied as a vertical force, we examine the sensitivity to source depth running simulations applying the source at greater depth (3613 m) in the layered velocity model (Figure 6).Comparing Figure 2b and Figure 6a we observe a significant reduction in path effects.To compare signal amplitudes with the shallow source, we plot in Figure 6b seismograms for receivers at position 1 and 8 from the left (compare first seismogram of Figure 6b with first row of Figure 3b and second seismogram of Figure 6b with first row of Figure 5b) showing a reduction in signal amplitude.Examining the sensitivity kernels, for instance, for station at position -2000 (Figure 6c), where a long wave train is observed, we can see how relatively weak later arrivals, even for the deep source, get trapped in the low velocity layer.However the effect is not as pronounced as for shallow sources that occur in low velocity layers.
Summarizing, we note that layers have a significant influence on the recorded traces producing a larger amplitude and an extended wave train in the seismograms, particularly for the shallow source.We note that the increase in signal amplitude is more pronounced for the horizontal component.We observe that the general trend is of increasing signal duration with increasing distance from the source.Furthermore, we see by means of the sensitivity kernels that, for the homogeneous case, the first seconds of seismograms are mainly controlled by the area between the source and the receiver in the shape of banana-doughnut (Tromp et al. (2005)), distorted by topography effects.Kernels plotted for the layered case show that the sensitivity area remains within the low-velocity layers.Although this is significantly more pronounced for the case where the shallow seismic source is located within the low velocity zone, it is also true for later arrivals for the deeper source although the seismogram coda are relatively weaker than for the shallow source case.The sensitivity kernel calculations demonstrate that waves stay trapped in the low velocity layers before extending over the entire surface of the edifice.

Sensitivity to source mechanism
In order to explore the sensitivity to the source type we replace the vertical force source with an explosive source.We show results of simulations considering the shallow source (113 m depth) Ricker wavelet with dominant frequency of 0.7 Hz. Figure 7 denotes the individually normalized velocity seismograms corresponding to the homogeneous model (Figure 7a) and the layered model (Figure 7b).We can see that the relative amplitude of the coda of the signals increase somewhat (compare Figure 7b with Figure 2b).Such observations are likely due to relatively weaker direct arrivals owing to the lack of direct S-waves in an explosive source.For the signal corresponding to the station at relative horizontal position 8000 metres (last station on the right in Figure 1a) we plot, in Figure 8, the sensitivity kernels finding the same result as before where trapped waves near the surface produce long coda, purely related to path effects.Qualitatively we see the same effects, independent of the source type.In the previous sections we investigated the effect of realistic near surface structure on seismic signals with central frequencies of 0.7 Hz, as it represents a typical value for LP seismicity.The effects are broadly similar for both a vertical force and an explosive source, therefore, we use only a vertical force as source mechanism in the following.

Conclusions
We examined the effect of the superficial layers on wave propagation in realistic volcanic setting.We performed systematic numerical simulations of wave propagation in a model We observe that topography and near surface low-velocity structure have a significant influence on the recorded traces as waves are trapped in near surface low-velocity layers.
These layers induce an extended wave train in the seismograms with a general trend of increasing signal duration with increasing distance from the source.In addition to distortion of the source signature we observe a significant increase in the amplitudes of horizontal component seismograms.This behavior does not depend on the source mechanism and is particularly acute for very shallow sources as seen on real LP data from Mount Etna (Bean et al. (2014)).Specifically, we observe strong tuning effects where 0.7 Hz and 1.0 Hz central frequency sources yield significantly different seismograms.Furthermore, results demonstrate that the seismogram coda depends on the properties of the whole near surface layer and not just on the portion below the station or along the source station direct path.
Results show the importance of the shallow velocity structure in shaping the recorded seismograms, implying that we need to account for the heterogeneities of the upper volcanic edifice when inverting for source mechanisms of shallow sources.Improving our knowledge of the near-surface velocity structures is thus a necessary step toward a detailed determination of shallow volcano source characteristics.

(
http://www.cgiar-csi.org/).The seismic network in the simulations consists of 15 recording stations distributed along the surface with a station spacing of 1000 m.Both shallow and deep sources are investigated.Two different velocity models have been considered.The first is a homogeneous one, with density ρ = 2384 kg m −3 , P-wave velocity v p = 3500 m s −1 and S-wave velocity v s = 2000 m s −1 .The second model is a heterogeneous one, and is obtained by combining the S-wave velocity structures at Mount Etna from Cauchie and Saccorotti (2013) for the top 360 meters and from Chiarabba et al. (2000) for the deeper part of the model.In the absence of clear constraints on Poisson's ratio (e.g.

Figure 1 :
Figure 1: (a) 2D model considered for the simulations.Recording stations are denoted in by black triangles and source positions by asterisks.Colours represent P-wave velocity structure for the heterogeneous case which is given in (b).(b) Seismic velocities and density structure for the heterogeneous model.The S-wave velocities for the first 360 m depth are from Cauchie and Saccorotti (2013) reporting the S-wave velocity structure derived from probabilistic inversion of Rayleigh-wave dispersion data.Deeper velocities come fromChiarabba et al. (2000).Corresponding P-wave velocities are obtained using a Poisson's ratio of 0.25 and densities using the Gardner's relationship(Gardner et al. (1974)).(c) Source time function (top row) and its amplitude spectrum (bottom row) used in the simulations.

Figure 2 :
Figure 2: Individually normalized velocity seismograms (horizontal and vertical component have been jointly normalized) corresponding to simulations in (a) the homogeneous, and (b) the heterogeneous model.The source is a Ricker wavelet with dominant frequency 0.7 Hz applied as a vertical force at 113 m depth.Synthetic seismograms are recorded at the black triangles in Figure

Figure 3 :
Figure 3: Waveforms and their spectral content corresponding to the horizontal component (left column) and to the vertical component (right column) of the synthetic signal recorded at relative position -6000 (first seismogram from the left in Figure 2) generated for a Ricker wavelet (f=0.7 Hz) applied as a shallow vertical force, using the homogeneous model (a) and the layered model (b)(See also Figure2).We observe that, for the homogeneous model the mean vertical amplitude of the signal is about 2 times the mean horizontal amplitude and that, while the vertical component increases by a factor of 2-3 times if a near surface layering is considered, the horizontal component increases by a factor of 7-8 times.

0. 7 Figure 5 :
Figure 5: Same as Figure 3 for the signal recorded at relative position 1000 (eighth from the left).

Figure 6 :
Figure 6: Simulations in the heterogeneous model.The source is a Ricker wavelet with dominant frequency 0.7 Hz (see Figure 1c), applied as a vertical force at 3613 m depth.(a) Individually normalized velocity seismograms (horizontal and vertical component have been jointly normalized).(b) Seismograms recorded on the stations located at relative horizontal positions -6000 (first seismogram) and 1000 (second seismogram).(c) Seismogram (top) recorded on the station located at relative horizontal position -2000 m (fifth station from the left) and P-wave (first column) and S-wave (second column) sensitivity kernels for each 3-seconds time slice.

Figure 7 :
Figure 7: Velocity seismograms (horizontal and vertical component have been jointly normalized) corresponding to simulations in the homogeneous model (a) and simulations in the layered model (b).The source is a Ricker wavelet with dominant frequency 0.7 Hz applied as an explosion at 113 m depth.

Figure 8 :
Figure 8: P-wave (left column) and S-wave (right column) sensitivity kernels for each 3-seconds time slice of the seismogram (top) generated by the shallow explosive source (0.7 Hz central frequency) in the layered model and recorded on the station located at relative horizontal position 8000 m.

4. 2 Figure 9 :
Figure 2b (generated by the vertical force with a 0.7 Hz central frequency) we can observe for the heterogeneous case that, despite the big spectral overlap (see Figure 1c) between both sources, simulations produce very different seismic signatures showing that the details of coda generation in the LP frequency range depend on the model structure relative to the seismic wavelength.
representative of a profile across Mount Etna, Italy, considering a homogeneous velocity model and a model with low-velocity layers at the top, as well as different source mechanisms and locations.We calculated the synthetic signals recorded by the stations along the surface and the corresponding P-and S-wave travel time sensitivity kernels by means the software SPECFEM2D.