Long Period (LP) Events on Mt. Etna volcano (Italy): the influence of velocity structures on moment tensor inversion

. We obtain very different solutions among the three models in terms of source location and mechanism decomposition. The overall shape of the retrieved source time functions are similar, but some amplitude differences arise, especially for the homogeneous model. Our work highlights the importance of including the unconsolidated surface materials in the computation of Green’s functions especially when dealing with shallow sources.


I N T RO D U C T I O N
The understanding of the origin of seismic signals on volcanoes is of fundamental importance to enhance our knowledge of volcanic systems and to monitor their activity. Volcanoes can exhibit a wide variety of seismic signal types (e.g. Chouet et al. 2006 and references therehein). Here we focus on long period (LP) seismic signals. LP events are characterized by low frequency waveforms  and are thought to have magmatic or hydrothermal origin (Chouet 2003). They are often considered to be associated with resonance of fluid filled cavities (Aki et al. 1977;Chouet 1996Chouet , 2003Neuberg & Pointer 2000;Nakano et al. 2003;Jousset et al. 2004) and their understanding is crucial to illustrate the shallow plumbing system of volcanoes. Recently, Bean et al. (2013) extended the observations of Harrington & Brodsky (2007) and proposed an alternative model for explaining shallow LP seismicity. They analysed the pulse-like nature of some LP events recorded on volcanoes and explained their origin as a slow failure of the weak shallow volcanic edifice close to the brittle-ductile transition. Their conclusions suggest that careful attention should be paid in using LP events as direct indicators of magmatic/hydortermal fluids.
An important tool to describe LP sources is moment tensor (MT) inversion (Kumagai et al. 2002Nakano et al. 2003;Jousset et al. 2004Jousset et al. , 2013Nakano & Kumagai 2005;Lokmer et al. 2007;Davi et al. 2010;De Barros et al. 2011). Many MT inversions on LP volcanic signals infer a tensile crack source mechanism (Kumagai et al. 2002Nakano et al. 2003;Eyre et al. 2013;Jousset et al. 2013). The orientations of these cracks span a full range from subhorizontal to subvertical. A different source mechanism was obtained by Davi et al. (2010) on Arenal volcano (Costa Rica) where their analysis of a LP signal recorded for an explosion led to an isotropic source mechanism. Although the inversion process itself is well established, many questions on uncertainties arise due to the lack of knowledge of the properties of materials traversed by seismic waves (Jousset et al. 2004;Bean et al. 2008). Although previously thought that kilometres long LP wavelengths 786 C. Trovato et al. are mostly insensitive to small variations of the volcanic structure, it has been shown that hundred meters thick superficial layers and features (volcanics above the sediments, dykes, magma pathways, etc.) can significantly degrade source inversion efforts (Neuberg & Pointer 2000;Kumagai et al. 2005;Bean et al. 2008;. This is because the complex stratigraphy of volcanoes has a strong impact on the seismic wavefield (Neuberg & Pointer 2000;Bean et al. 2008). Nevertheless, many inversions were conducted in a homogeneous half-space (Ohminato et al. 1998;Legrand et al. 2000;Kumagai et al. 2002;Jousset et al. 2004;Lokmer et al. 2007;De Barros et al. 2011). There were also efforts to invert the source using a more complex media, such as a two-layered medium  or a heterogeneous medium (Davi et al. 2010;Jousset et al. 2013;Eyre et al. 2015); however, little is mentioned about the quantification of the introduced error in the case where the model is completely or partially incorrect. At present, detailed velocity structures of volcanoes are still rarely available (Davi et al. 2010). Here we extend the work of Bean et al. (2008) on the influence of the velocity structures on moment tensor inversion. We analyse the effect of different geological models with increasing complexities regarding their influence on successful source inversions.
We focus our attention on one of the most studied volcanoes in the world, Mt Etna, located in eastern Sicily island (Italy), the largest active volcano in Europe. More than 600 000 people live nearby this active volcano (Chiarabba et al. 2000). It covers an area of about 1250 km 2 and reaches a maximum elevation of ≈3330 m. It's characterized by almost continuous eruptive activity from its summit craters and frequent lava flows from fissures opened on the flanks . Due to the high volcanic activity many seismic signals are now recorded continuously . Since 2003 the permanent network of broadband stations has been installed by the INGV (Istituto Nazionale di Geofisica e di Vulcanologia) and LP events have been addressed in a number of studies Saccorotti et al. 2007;De Barros et al. 2009, 2011Cannata et al. 2009Cannata et al. , 2013. They appeared in periods of quiescence or during unrest episodes Saccorotti et al. 2007). They are often difficult to be distinguished from the sustained volcanic tremor accompanying eruptions ), but they may not be directly related to the eruption processes De Barros et al. 2011). The mechanism of the LP events suggested resonating phenomena at a relative shallow depth (∼300-1200 m. below the summit) De Barros et al. 2011). Recently, Bean et al. (2013) proposed a new model for explaining the shallow LP seismicity recorded in occasion of the 2008-2009 eruption of Mt Etna (De Barros et al. 2009, 2011. They recognized that, while summit stations recorded pulse-like low-frequency signals, the same records on further stations appeared as classical resonating LP signals. They attributed the apparent resonance of these low-frequency seismic events to propagation effects and not being source related. Their model hypothesizes that those LP events are consequence of failure in materials close to the brittle-ductile transition. The brittle-ductile transition in shallow volcanic material is not supposed to be related to high temperature and pressure, but to the low friction angles of the unconsolidated shallow volcanic deposits. Similar conclusions have been drawn by Eyre et al. (2015) for Turrialba volcano, Costa Rica.
In this study, we carry out a synthetic inversion verification test ('blind test'). We build four different structure models with increasing geological complexity. We suppose that the fourth most complex model corresponds to reality (state of reality) and that the three other models correspond to the best knowledge we have of the geological properties of the volcano (state of knowledge). Hence, we compute Green's functions for the first three models and synthetics data in the fourth most complex model for a tensile crack source mechanism. Thus, we will perform three moment tensor inversions of this realistic synthetic dataset, using the Green's functions database outlined above. We vary the number of free parameters in our inversions, as well as the number of receivers. For the given synthetic scenario, we will discuss how well the inversion procedure can reproduce the original location, mechanism and source time functions (STFs) using the different structure models and stations configuration. The use of synthetic scenarios allows us to extend the previous studies and to systematically assess the influence of the velocity model complexities to the accuracy of the retrieved source mechanism. We then extend the work of De Barros et al. (2011) analysing a real LP event by means of the three geological models with variable complexity.

Velocity models
Geological mappings of Mt Etna volcano have been performed since more than a century. Surface units have been mapped by De Beaumont (1836) and the first geological maps of Etna volcano were published in the 19th century (Waltershausen 1844(Waltershausen , 1880. In the last decades, official geological maps were updated twice (Romano et al. 1979;Branca et al. 2011) and many geological surveys have been carried out to map deposits along the steeps of Valle del Bove (Calvari et al. 1994;Coltelli et al. 1994) and integrate in the Italian geological map of the surroundings (Pasquare et al. 1992;Branca et al. 2009Branca et al. , 2011. On the other hand, many geophysical seismic surveys have been also carried out (Hirn et al. 1991;Cardaci et al. 1993;Luca et al. 1997;Villaseñor et al. 1998;Chiarabba et al. 2000;Laigle et al. 2000;Patane et al. 2002;Cristiano et al. 2010;Cauchie & Saccorotti 2013) analyzing the velocity structure properties of the edifice. Following this studies we prepare four different models of the velocity structure by increasing its complexity on the depth variation. All the models constructed with help of a meshing tool (CUBIT-13.2 from Sandia Laboratories) include topography from the Digital Elevation Model (DEM) of Mt Etna with a 50 m spacing. Horizontally we prepare a model extending 19.6 km in the EW and 16 km in the NS direction, with a max height of about 3300 m (Fig. 1), large enough to minimize reflections from the model boundaries.
We use four (4) different models shown in Fig. 2. Model (S1) is homogeneous with P-wave velocity of 2000 m s −1 and Vp/Vs ratio of 1.73 (the value of 2000 m s −1 is taken from De Barros et al. 2009). The second model (S2) takes into account a low-velocity surface layer of 300 m thickness inferred from Bean et al. (2008). The third (S3) and the fourth (S4) models are more complex. In these two models we adopt the gradient model of Mt Etna according to the geological map of Branca et al. (2011). The strong stratigraphicgradient of the volcano is represented by different piled layers, thus are characterized by topography shape and become flatter with depth (towards a proportional smoothing function) until the sea level. We define the velocities at depth according to Chiarabba et al. (2000), Branca et al. (2009), Cristiano et al. (2010 and Patane et al. (2002). Model S3 has a homogeneous surface layer of 300 m (Fig. 2). Model S4 (Fig. 2) has a strong gradient structure in the shallow depths down to 360 m as inferred from Cauchie & Saccorotti (2013). P-wave velocities (Vp) are derived from S-wave velocities (Vs) as, V p = √ 3V s assuming a constant Poisson's ratio of 0.25. The material densities are computed following the formula proposed by Potter (1998) in function of Vp. We use models S1, S2 and S3 to compute Green's functions for the inversion while model S4 is used to compute the synthetic data.
Before interpreting the full waveforms computed in 3-D complex geological media including sharp topography and non-planar layering, we look at the response of 1-D layered structures consisting of horizontal layers overlying a homogeneous half-space with same profile as in Fig. 2. In addition to these geometric simplifications, the details of the seismic source are avoided by considering a plane SH-wave propagating vertically from the half-space (i.e. no incidence angle has been considered). Under these hypotheses, only reflection/transmission of the SH waves will occur at the interfaces; however, the interpretation of such simple configurations will give us valuable information about the impact of the 1-D soil layering on vertically propagating SH-waves. For solving the 1-D wave propagation in our four models S1 to S4, we compute the semi-analytical solutions in the frequency domain by the Thomson-Haskell propagator matrix method (Thomson 1950;Haskell 1953) and use the inverse Fourier transform to obtain the time domain waveforms. The plane wave is injected at 3 km depth and the impulse STF has a flat frequency spectrum up to 10 Hz. Fig. 3(a) shows the filtered (0.2-2 Hz) velocity waveform at the free surface. The waveform S1 has the exact shape of the STF multiplied by a factor two. The original shape of the STF is almost unchanged for models S3 while models S2 and S4 are significantly subjected to reflections/transmission effects due to the higher velocity contrast of the shallow layers. Model S3 shows arrival times comparable to the reference model S4, while those in models S1 and especially model S2 are considerably delayed (∼1-1.5 s). Recorded amplitudes are higher for model S4 decreasing towards models S1. The frequency content of velocity traces (Fig. 3b) shows a single peak for model S4 (f ∼ 1.2 Hz) while the other models show different peaks with the main energy focused at lower frequencies. This simple comparison is obviously not representative of the complex Mt Etna geological context with the presence of topography, but it clearly illustrates that models S1, S2 and S3 are different enough from model S4 to justify the proceeding of the synthetic test.

Methodology
In an elastic medium, the n-th component of the displacement (u n ) at a point x at a time t is given by the convolution between the source-time function of the moment/single force and the medium response, that is Green's functions (GFs; modified from Aki & Richards 2002): where M pq is the pq-component of the seismic moment, F p is the single force acting in the p direction and G np represents the medium response (Green's function) for the nth-component displacement due to a unit single force F p and G np,q means the spatial derivative with respect to the q-component at the source position. The asterisks indicate an operation of convolution and the Einstein summation convention is applied. The numerical method used to compute the synthetic seismograms is the spectral-element method (SEM). This method has been developed in computational fluid dynamics by Patera and Maday (Patera 1984;Maday & Patera 1989). It was then introduced in computational seismology a decade later to compute with accuracy wave propagation in complex geological media (Faccioli et al. 1997;Komatitsch & Vilotte 1998). The SEM is a high-order finiteelement approximation in which the consistent choice of an orthogonal polynomial basis and of a Gauss numerical quadrature leads to the convergence properties of spectral methods. Because of the use of high-order piecewise polynomials basis, numerical dispersion is significantly reduced compared with the classical finite-element method (FEM;De Basabe & Sen 2007;Seriani & Oliveira 2008). The reader is referred to Komatitsch et al. (2005) and Chaljub et al. (2007) for review articles presenting the numerous developments of SEM, and to Moczo et al. (2014) for a historical presentation and recent applications to seismic wave propagation in alluvial valleys.
The FEM, and therefore the SEM, is particularly well adapted to compute wave propagation in media where the relief features (such as mountains, hills, creeks or volcanoes) are present because the free surface condition (and more generally the continuity of traction) is said to be a 'natural boundary condition'. In other words, the free surface, no matter what shape it has, is accounted for in the weak form of the equations to be solved and does not have to be explicitly enforced at the elemental level. This allows surface topography to be accounted for in SEM, as long as the elements can honor the shape of the free surface without any aliasing. For this study, a special care has been taken for the generation of the meshes so that the finite elements of size 50 m at the free surface follow the volcano's topography provided by a digital elevation model (DEM) of resolution 25 m.
In this article, to calculate GFs in the elastic medium with irregular surface topography and synthetic seismograms, we use the open-source code EFISPEC3D (http://efispec.free.fr). This computer program (under double licenses CeCILL-V2 and GNU-GPL-V3) solves the 3-D equations of motion using a continuous Galerkin spectral-element method. The correctness of the implementation of the spectral-element method into this code has been thoroughly verified in De Martin (2011) and Chaljub et al. (2015). EFISPEC3D is used in computational seismology to better understand the impact of lithological and topographic effects on near-surface Green's functions (Maufroy et al. 2015).
We put potential source positions within a volume of 1.000 × 1.000 × 800 m 3 located below the main crater of Etna volcano between 2.2 and 3 km a.s.l. (Fig. 1). Among the 28 receivers used in this study, 13 receiver locations correspond to the stations of the permanent network operated by INGV, 12 are from temporary surveys (Saccorotti et al. 2004;Lokmer et al. 2007;De Barros et al. 2009) and 3 receivers (synthetic stations) are added to guarantee the azimuthal coverage for our synthetic test (Fig. 1). In order to treat a large number of source locations we take advantage of the reciprocity (Aki & Richards 2002) to calculate GFs.
We carry out the inversion in frequency domain for eq. (1), which is schematically written as a vector-matrix equation: where u is the data matrix, G contains the Green's functions and m is the moment tensor and single forces components that we aim to obtain. We perform the inversion for the model parameters m without applying any a priori constraints to the solution (hereafter called 'unconstrained inversion'). We define the misfit (R) function as:

790
C. Trovato et al. where superscript T denotes a transposed matrix. The least-squares solution of eq. (2) is given by (e.g. Menke 1989): This inverse problem (eq. 2) can be solved either for six independent moment tensor components (MT) (assuming no single forces), or six moment tensor plus three single forces (MT + F). The inversion is carried out for each position of the source (14 196 positions at 40 m spacing). Comparing the value of the misfit R from each inversion we can estimate the best fitting source position. For analyzing the estimated solution in terms of their mechanism, we use the principal component analysis (PCA) through a singular value decomposition (SVD ;Vasco 1989;. This technique assumes the existence of a unique STF for all the six components of the moment tensor (see Vasco 1989 for further details). We then decompose the moment tensor solution into isotropic (M ISO ) and deviatoric (M CLVD + M DC ) parts after Jost & Herrmann (1989) and Vavryčuk (2001).
Additionally, we also perform a constrained inversion following the approach by Lokmer et al. (2007), assuming either a tensile crack or isotropic source mechanism. Eq. (2) is rewritten as where f is a function of strike φ measured from the North in the clockwise direction and dip θ, independent of frequency. Our inversion reduces to finding a single parameter M 0 (ω). We perform the grid search spanning from 0 • to 360 • for strike (φ) and from 0 • to 90 • for dip (θ ), every 10 • for the tensile crack mechanism. For an isotropic source one inversion is enough as the function f has a unique expression. The inversion procedure has been verified for different source mechanisms by producing synthetic data in the velocity model S2 and inverting with Green's function for the same velocity model. MT inversion and SVD delivered the perfect solution validating our implementation.

Synthetic data
Previous inversions of LP signals on Etna volcano De Barros et al. 2011) suggest quasi-vertical crack orientations. Hence we simulate, as the synthetic source mechanism, a point source of a vertical tensile crack (φ = 45 • , θ = 90 • ) at two different depths located below the summit craters: at 2880 m.a.s.l. (shallow source, ∼400 m depth) and at 2240 m.a.s.l (deep source, ∼1.2 km depth). We use a Ricker wavelet as STF with the main energy in the frequency range 0.2-1.2 Hz (typical for LP events, Chouet 2003) and an amplitude of 4 × 10 10 Nm. As already mentioned, model S4 is used to calculate the data. Time step is t = 1 × 10 −4 s, for a duration t tot = 20 s. A single simulation of 2.3 × 10 6 hexahedron elements takes about 18 hr on 192 cores on our local server (AMD Abu Dhabi at 1.6 GHz). Synthetic data are computed without adding noise.
So, for each velocity model we calculate three inversions, MT (moments only), MT + F (moments plus single forces) and CON-STR (constrained inversion).

R E S U LT S
We first carry out the unconstrained inversion in order to investigate the reliability of the solution and the uncertainty between the different velocity models.
In the following, we will discuss the source location, the source mechanism and the STF obtained from MT/MT + F inversions with models S1, S2 and S3. The 12 stations located nearby the summit, offering a proper azimuthal coverage, are used (Appendices A and B). We choose to include only 12 stations located in the nearintermediate field, as we want to mimic realistic number of available stations.

Source location
First, for the two given source depths (∼400 m and ∼1.2 km depth), we explore how well the inversion procedure can retrieve the true position in different structural models. We evaluate the misfit from the moment tensor plus single force (MT + F) inversions. The comparisons of the best fit waveforms for the shallow and deep sources are shown in Appendices A and B, respectively. The overall shape of the original signals is well reproduced in all three velocity models. The stations closer to the source reproduce better the original STF shape as, for small distances, the wavefield is less subjected to attenuation, scattering and reflection phenomena. This is the case for stations et08, et06 and et09, while differences on the waveforms are stronger for farther stations (e.g. cl01 and et99).  Tables 1  and 2. For the shallow source, the minimal misfit found in model S2 coincides with the original source position. For models S1 and S3 the obtained source locations are shallower (∼100 m) than the original one, at the upper limit of our parameter search. The synthetic test has been built as a 'blind test', and a further computation of Green's functions is computationally expensive. Hence, we consider these locations as a local minima and not a global one.
For the deep source we get good (∼200 m distance from the original position) horizontal and vertical resolution for both S1 and S2 models. For model S2 the location is slightly better constrained, as the lowest misfit value sharply converges to a single position. For model S1 we can observe a large spreading of low misfit values, even if the lowest misfit still points to the right source position. The reason is that the shape of the waveform in a homogeneous model is very weakly affected by the perturbations of the deep source location around its true position. Model S3 points to a quite distant (∼400 m from the original position) source location and the value of misfit (Table 1) is considerably higher than the two other models.
Epicentral locations both for the shallow and deep sources are very well constrained except for the deep source in the velocity model S3. The good radial/azimuthally station coverage guarantees a good resolution in retrieving the original epicentral position. The non-negligible influence of the strong velocity gradient in the shallow part of model S4 is better resolved when considering simple velocity models S1 and S2. On the opposite, especially for the deeper source, the velocity contrast reproduced in model S3 degrades the solution. This implies that a simple model is better to use if we do not know exactly the velocity structure. This is particularly evident for the deeper source where lateral heterogeneities play a major role (Appendix B). On the other hand, vertical locations are less resolved and appear more sensitive to the wrong velocities definition.
In summary, source location towards MT-Inversion (and especially vertical resolution) is very sensitive to the choice of the velocity model. The results obtained here, even for such a simple case (there is no noise and the reference model S4 is still simple if compared to reality) point out that other locations methods such as   (left-hand panel) and deep (right-hand panel) source obtained in the inversion of minimum R using different structure models. Axes represent the relative source location with a 40 m spacing. The real source position is represented by (0,0,0). From top to bottom velocity models S1, S2, and S3 used for the inversion. Lowest misfit value is represented by the synthetic slices intersection for each velocity model. amplitude decay, cross-correlation coefficient and semblance (Cannata et al. 2013 and references therein) less sensitive to velocity definition should instead be used.

Source mechanism
For the best source position obtained above, the source orientation and isotropic/deviatoric decomposition for each model are estimated in Table 1 (shallow source) and Table 2 (deep source). We compute as well the moment magnitude of the selected event as: where M 0 is the highest absolute retrieved seismic moment after SVD. Table 1. Summary of the results in MT + F and MT inversions for the shallow vertical crack, using different structure models (S1, S2 and S3), respectively. For the best misfit in each inversion, the fault mechanism (strike, dip) are calculated and the decompositions are performed. We also report the validation misfit between the observed MT solution and the original one used to reproduce the tensile crack in the synthetic test.

Shallow source
The crack strike (φ) is well retrieved for both inversions (MT + F and MT) while the obtained crack dip (θ ) is close to a solution of a horizontal crack rather than a vertical one for the solution including single forces. The minimum misfit R is found for model S1. In terms of isotropic/deviatoric decomposition, the MT + F inversion points to the right ratio, letting the predominant component in C ISO , C DC being close to zero and C CLVD showing values close to the given one for all three models. The MT inversion without single forces points to a very low C CLVD value and a high C DC . In this case the mechanism would be interpreted differently with a strong double-couple component and a mixed ISO/DC solution. Table 2 shows the results for the deep vertical crack. Again the lowest misfit values are found for S1 model (both MT + F and MT solutions). Model S3 shows a very high misfit value. The crack strike (φ) and dip (θ ) angles are very similar to the true ones for all the three models and the best solutions are obtained for model S1 (MT + F) and S2 (MT + F). In terms of mechanism, the decompositions for models S1 and S2 give a low C DC component and the crack solution is well retrieved, while model S3 tends to overestimate the C DC component (16 per cent for the MT + F solution). The same observations are brought for the MT inversion except that we find higher C DC component contributions.

Deep source
In summary, we find that model S1 gives the lowest misfit value in both the MT + F and MT inversions for both source depths. The crack orientation is better retrieved by model S2 such as the isotropic/deviatoric decomposition. For the shallow source it is difficult to retrieve the right crack dip angle with any of the velocity models in the MT + F solutions; the shallow vertical crack might be interpreted as a shallowly dipping one. The dip angles are better retrieved when considering MT solutions, but we are still far from the right orientation angle. Moreover, when not considering single forces, we find very high C DC values which complicate the interpretation of the solution. The magnitude of the event is quite well retrieved by all solutions (Table 1): slightly overestimated for MT + F (especially model S1) and slightly underestimated for MT-only (except model S1). The strong retrieved single forces have the effect to increase the strength of the seismic moment. This is probably due to interference of the radiated seismic waves by both the MT and single forces On the opposite, the deep source orientation is well retrieved by all the three velocity models for both solutions MT + F and MT-only. The retrieved values for the magnitude of the event (Table 2) are equivalent between MT + F and MT-only solutions (slightly underestimated in respect to the expected event magnitude) and the appearance of spurious single forces does not influence the energy of the radiated seismic waves. In first approximation, the 'blind test' tells us that little matters the chosen velocity model, we encounter problems in defining the right source orientation for shallow sources, but for deep ones things turn out to be better and we always obtain a good solution. The inability in retrieving the correct solution for the shallow source could be due to the strong velocity contrast just below the ground surface of velocity model S4. The shallow crack source is embedded inside these low velocity sector, hence none of the three Green's functions models is able to model correctly the strong impedance contrast of model S4 (especially near-field terms, Lokmer & Bean 2010) and this eventually leads to errors which tend to be accommodated by single forces which, in turn, degrade the MT solution. The same does not apply for the deeper source because of the longer paths followed by waves and predominance of intermediate-and far-field effects.
In order to explore the influence of the strong velocity contrast of model S4 on the retrieved MT components for each velocity model, in the next section we compare our retrieved STFs for each MT component with the expected ones.
where MT R is the retrieved MT solution, MT T is the original MT tensor and N is the number of time series.

MT + F
In the MT + F inversion the force terms show high amplitudes especially for model S1 (Fig. 5) (Tables 1 and 2) between the original and the retrieved STFs (considering only the MT terms) are very high for the three models. Model S1 shows the highest validation misfit (28.972). Model S2 shows the best match with the original solution. The original STF is generally well reproduced (especially models S2 and S3), but the M zz component suffers of overestimation in amplitude thus leading to the high validation misfit and the wrong source mechanism orientation interpretation. For the deep source, the amplitudes of the seismic moment are always underestimated. The validation misfits are considerably lower than for the shallow source with model S1 showing better correspondence with the original STF. Anyway, the overall shape of the retrieved STF is similar and well reproduced for all different models.

MT-only
The MT inversion (Fig. 6) shows similar results, but the validation misfits are considerably lower than the solution including forces and the original STFs are better reproduced. Again, the Mzz component of the shallow source does not match the true solution especially for the simplest model S1. This is due to the quasi-horizontal layering in which the wave conversions occur. Models S2 and S3 show comparable validation misfits (∼0.85) and the amplitudes of the moment tensor components are comparable to the original STFs. For the deep source the misfit values are lower, the same as for the MT + F solution. Best matching between the original and the retrieved solution is obtained by model S1 (validation misfit = 0.545). The overall amplitude of the STF is in general underestimated, likely because the velocity at depth is much lower in S1 than in the true model S4.
In summary, the synthetic test shows that the deep source mechanism is correctly retrieved by both the MT + F and MT solutions. The shallow source, on the other hand, suffers from large errors in the retrieved STF which strongly influence the mechanism decomposition (as highlighted previously in . Main uncertainties for the shallow source are due to the large leakage between the true and retrieved M zz component. This implies that we are unable to correctly resolve the isotropic component of the MT. Consequently, the MT + F solutions and the MT solution for model S1 would be interpreted as a shallowly dipping tensile crack. The MT solution for models S2 and S3 point to a quasi-vertical crack as expected, but the appearance of the non-existing double-couple components complicates the mechanism interpretation.

Constrained inversion
Finally we perform the constrained inversions. Our input source is a vertical tensile crack oriented 45 • clockwise from North. We perform the inversion assuming a tensile crack and an isotropic source, respectively. The results are summarized in Table 3 and the Misfit results for the crack mechanism are shown in Appendices C and D. All models show a good solution with the lowest misfit indicating a tensile crack mechanism and pointing to the correct angles [strike (φ) and dip (θ)] orientation for both the shallow and the deep source in the given parameter ranges. Figs 7 and 8 show the comparison between the original and retrieved STFs for all three models and depths from the inversions with and without single forces, respectively. The retrieved amplitudes in models S1 are overestimated in all the inversions. The amplitude is twice than expected for the shallow source. STF shape is well retrieved in both models S2 and S3 for the solutions with and without singles forces. For the deep source a phase shift between the original and retrieved STF occurs due to travel time delays caused by different velocity model. Generally both models S2 and S3 offer a good solution in both angle pairs and STF, while model S1 tends to overestimate the STF amplitude.
We now apply our results obtained by the synthetic test to the inversion of a real event recorded at Mt Etna.

R E A L C A S E : A N L P E V E N T I N 2 0 0 8
Despite our synthetic test has been designed to mimic reality, in the real world MT inversion is subject to significant uncertainties which should strongly influence our ability in retrieving the correct solution. Here we want to show the influence of the choice of a particular velocity model on the inversion process, thus we carry out an inversion of an LP event recorded on Etna in 2008 during a high resolution seismic survey (De Barros et al. 2009, 2011. The considered event was recorded on June 19, 2008 and belonged to the second family of events identified in De Barros et al. (2009). The source mechanism was analysed (De Barros et al. 2009Barros et al. , 2011 (i) by locating the event with a time delay technique based on cross correlation and (ii) by identifying the source mechanism performing a MT inversion using a homogeneous model (same as our model S1). The mechanism was retrieved as a subvertical crack oriented φ = N340 • E and inclined θ = 50 • (see De Barros et al. (2011) for further details). While De Barros et al. (2011) used 16 stations in their inversion, we choose 12 stations with a good azimuthal distribution around the source (Appendix E) in order to reproduce a context similar to the one chosen for our synthetic test.

Unconstrained inversion
Appendix E shows the comparison between the original filtered data and our synthetics resulting from the MT + F inversion for the three models separately. Stations etsm and et08 show the highest amplitude signal and thus contribute more to the final solution. For these two stations we get a good correspondence between the observed and the retrieved signals for all three models. Farther stations do   Barros et al. (2011). The MT + F solution from models S1 and S3 also varies, in particular in terms of strike. The dip of θ = 67 • from model S3 is comparable to the one found by De Barros et al. (2011). Fig. 10 shows the retrieved STF for the three models after the MT inversion. Here the MT solution in model S1 shows a higher amplitude than the two other models. As MT and MT + F solutions are often used to estimate the volume of fluids or gas mobilized at the source (Ohminato et al. 1998;Hidayat et al. 2002;Davi et al. 2010;Jousset et al. 2013), it is obvious that the results may be very uncertain, depending on the discrepancy of the used velocity model from the true one. The amplitude difference is also present in MT + F case but it is less remarkable. The overall shapes of the retrieved STFs for the three models are quite similar and model S2 shows the simplest solution. For all the three models, the diagonal of the moment tensor is largely dominant while non-diagonal elements show much lower amplitude. Non-diagonal elements in the solution without single forces tend to be overestimated compared to the solution including forces.

Constrained inversion
The high C ISO components retrieved for the unconstrained inversion for the three velocity models suggest a possible isotropic source mechanism and so does the constrained inversion which shows slightly lower misfit values for the explosion solution (MT + F, Table 5). It is worth noticing that a tensile crack with a low Poisson's ratio in the source region could also lead to the same eigenvalues ratio. In terms of the orientation of the crack solution, the parameter search for the MT + F solution does not indicate a clear orientation (Appendix F). The MT + F solution shows a narrow range of misfit values spanning from 0.44 to 0.5 for models S1 and S2, and from 0.5 to 0.6 for model S3. The minimum misfit found for model S2, MT + F solution, shows a strike orientation similar to the one found by De Barros et al. (2011), but the inclination of the fault once again results in a subhorizontal instead of a subvertical crack. The MT-only solution results more stable, and all three models point to the same solution (φ = ≈300 • and θ = ≈50 • ), but the misfit values are very high (≈0.8). As the misfit values for the constrained inversion fall in a very narrow range, the solution is subject to difficult interpretations, that is it is difficult to discriminate between an isotropic and a tensile crack source mechanism.

Q UA N T I F Y I N G T H E ' G O O D N E S S O F S O L U T I O N '
In this section, we quantify the reliability and sensitivity of the solution in our inversion framework. In order to achieve this task, the Green's functions should be repeatedly calculated for a large number of perturbed velocity models. However, in 3-D heterogeneous medium with topography this would be extremely computa-tionally expensive, hence we use a different approach instead: for each velocity model, we perform a large number of inversions for randomly chosen network configurations. In this way we implicitly include different parts of velocity models into inversions. From such a large number of solution for a particular velocity model we calculate (i) robustness of the inversion for a given model, and (ii) departure of the retrieved from the true model. In this test, all MT inversions are performed for the fixed source location at the original point in order to get rid of errors due to mislocation and tackle directly the influence of the velocity model on the source mechanism. We select 21 seismic stations located at distances shorter than 5 km from the source. We then randomly select 8-16 stations from the available 21. In this way we obtain 1350 station combinations and, for each of them, we perform MT-inversion for all the three velocity models and both source depths. After applying PCA to the MT solutions we analyse the source properties and orientation as outlined in Section 2.2. For each velocity model/source combination, we compute the median and the absolute median deviation. In order to make the results more intuitively comprehensible, we also compute the slip direction α after Vavryčuk (2001) representing the angle between the fault plane and the slip vector, that is α = 90 • for a pure tensile mechanism, while α = 0 • for a pure shear faulting.
The results of MT + F and MT-only inversions for both the sources are reported in Fig. 11. The results confirm many of the conclusions outlined during the previous synthetic test (Section 3). Model S1 always delivers the lowest misfit in both MT + F and MT-only inversions, but for the shallow source it also leads to the worst validation misfit in both the MT + F and MT-only inversions.
More generally, we find that the MT + F inversion (Figs 11a and c) is always able to retrieve the correct strike ( ) for both the shallow 798 C. Trovato et al. Figure 9. Misfit with respect to the source position for the LP event. The location through MT inversion is conducted in three different velocity models (S1, S2 and S3, respectively). The axes in meters represent the relative distance to the original source location determined by De Barros et al. (2011). Original solution (x, longitude, 49 950 km, y, latitude, 4 178 450 km, z, height, 3 km). and deep sources, but it fails to correctly determine the dip (θ) (θ ≈ 45 • ) for the shallow source and well retrieves it for the deeper one. Finally the angle α is well retrieved for the shallow source, but is highly underestimated for the deeper one (α ≈ 50 • ), that is the mechanism would be misinterpreted as a mixed tensile/shear source.
When we focus on the MT-only solutions (Figs 11b and d), both the strike and dip orientations are well retrieved for both the sources even if the dip angle for the shallow source is affected by high uncertainties. We also find that the angle α is underestimated for both the sources and the results show again high uncertainties (especially for the shallow source). The C DC component is very high (C DC ≈ 40 per cent) for both the sources.
In the following paragraph we compare the validation misfits (expressed as L1-norm departure between the retrieved and the true solution) for each MT component (Fig. 12a). Looking at the MT + F solutions (Fig. 12a), the shallow source shows the high validation misfits for all the components. The worst result is obtained in M zz component, especially for velocity models S1 and S2, and this leads to the wrong dip angle retrieved after MT decomposition. Model S3 shows lower validation misfits in all the MT components for the shallow source, but the MT decomposition still points to the wrong source orientation. The deep source shows lower validation misfits and the three models deliver comparable (and acceptable) results. Compared to the MT + F solutions, the validation misfits from the MT-only inversion (Fig. 12a) are considerably lower for both the shallow and the deep sources. In addition, the MT-only solution is less sensitive to the selected subset of stations (smaller fluctuations around the median solution). In this case the worst result is obtained again for the M zz component from model S1; also, the off-diagonal components show quite a large departure from their true values for the shallow source. Generally, including single forces in the solution degrades the match between the observed and the retrieved STFs. M zz is the most affected component and simple models (S1 and S2) deliver the worst results.
On the other hand, in order to figure out the stability of the inversion process, we plot (Fig. 12b) the validations misfit against the number of the receivers used in the MT inversion. The validation misfit tends to decrease when increasing the number of stations. This is particularly evident for the shallow source (both the MT + F and MT-only solutions) and for the homogeneous model S1 (The validation misfit is double when using 8 instead of 16 stations). However, we do not see any significant improvement in the retrieved source parameters (not showed here). The convergence of the results does not mean that we get the original source mechanism more closely by simply playing on the stations network. For the shallow source, including 16 stations (realistic case for real LP events) in the MT inversion is not enough to obtain acceptable results. In all the synthetic tests, we always misinterpret one of the source parameters, either dip or α angles. More generally, none of the velocity models leads to the correct solution (shallow source), while the three velocity models deliver very similar results.

Implications of the results
We performed synthetic tests to investigate the sensitivity of seismic source inversion results to the choice of the structural model. The homogeneous model (S1) shows lower misfit values than the other complex velocity models, but the retrieved STFs strongly deviate Downloaded from https://academic.oup.com/gji/article/207/2/785/2583626 by guest on 17 December 2021  Barros et al. (2011) is also shown at the bottom. MT-only inversion results for model S2 are missing as the retrieved STFfor different MT components were not represented by a unique source-time history (more than one significant singular values present in the SVD decomposition of the solution). This is normally due to inability of the MT-only solution to account for the uncertainties in the velocity model (it can be often observed when performing MT inversion of real data).  Fig. 2). Finally model S2 offers the best result in our synthetic test. It is worth nothing that the lowest misfit values obtained for the shallow source with model S1 correspond to the highest validation misfit values computed between the observed and retrieved STFs, leading to an important conclusion that the smallest inversion misfit does not always corresponds to the most accurate source mechanism (the same conclusion as in Bean et al. 2008). In the synthetic tests the results are similar and an approximately good solution can always be retrieved for each considered velocity model (i.e. the shallow source mechanism is well retrieved when a constrained inversion is performed; for the deep source we get a good solution even for an unconstrained inversion). Then we analysed a real event recorded on Mt Etna during a high resolution seismic campaign in 2008. Our results show that, in this particular case, the model with the lowest misfit value is the surface layer model S2. The analysis of a real LP event highlights the influence of the choice of a particular velocity model on the retrieved solution, that is the interpretation of the source mechanism varies for each considered velocity setting. These differences in convergence of the solutions between the synthetic and real data sets could arise because in the synthetic tests we are dealing with a simplified version of the reality. In particular, we ignored the following facts: any noise in the signal, complex source process different from a simple tensile vertical crack and further heterogeneities in the velocity model (especially lateral heterogeneities). All these factors are common for most MT inversions of LP events performed everywhere and could be source of errors in the retrieved MT solutions.

General remarks
The inversion tests were performed using the topography of the Etna volcano, but the overall remarks can be taken into account in any MT inversion for real LP events at volcanoes. The summary of our observations is given below: (1) Location: Locating events by the MT inversion grid search may lead to an incorrect source location and its mechanism (see results for model S3 in Section 3). Thus, we suggest to locate the events with some alternative technique, less sensitive to the choice of velocity model.
(2) Inclusion of single forces: De  suggest that the inclusion of single forces into MT inversions accommodates er-rors arising due to the mismodelling of the structural properties of the volcanic edifice. This is demonstrated by showing the equivalence of the radiation pattern of the vertical force in a homogeneous medium and that of the converted S waves at the low-velocity interface inside volcanic edifice. However, their examples apply for a flat medium without topography. Our results here show that although the inclusion of single forces can indeed accommodate the mismodelling in certain cases, it generally increases the discrepancy of the retrieved from the true solution. In addition, there is a significant energy leakage between the vertical force and the Mzz component, so we suggest to generally avoid the MT + F inversions.
(3) Shallow source: The solutions for the shallow source are strongly influenced by the generation of the surface waves and converted phases present in the waveforms calculated for model S4, hence the inversion is subjected to high uncertainty and misinterpretation. Adding complexities in the velocity model used to compute Green's functions does not necessarily lead to the correct solution. Generally better results are obtained without single forces, but this leads to spurious double-couple components. Even without the inclusion of single forces, the Mzz component for shallow sources shows increased sensitivity to the shallow part of velocity models.
(4) Deep source: Deep source is well retrieved by all the models for not constrained and constrained inversions.
(5) Constrained inversion: Good results for the orientation of the shallow crack can be obtained for every velocity model only when performing a constrained inversion (same conclusion as in Bean et al. 2008 andLokmer et al. 2007). On the other hand, when dealing with real data the constrained inversion may deliver ambiguous results. If a converged solution cannot be found when performing a constrained inversion, we should move to different strategies [inversion of tilt components (Maeda et al. 2011;Chouet & Dawson 2015;Thun et al. 2015;van Driel et al. 2015), lower frequencies, etc.].
(6) Velocity model: An important question addressed in this study is whether including complexities in the velocity model improves our ability in retrieving the original source mechanism. The synthetic tests show that, for the deep source, the three velocity models deliver similar and acceptable results. However, the model S1, with the underestimated velocity at depth, overestimates the STF. It suggests that in the case of deep sources, the best available tomographic model should be used. One the other hand, for shallow sources, the vertical component (M zz ) of the MT tends to be incorrectly retrieved by any velocity model. A likely reason for this is improperly captured surface waves and pronounced converted phases present in the wavefield. This makes the results for shallow sources unreliable and directly affects our efforts to estimate the amount of gas/magma    involved in the source processes. To summarize, for deep sources we are allowed to use a simple velocity model of the volcano (if a comprehensive description of the geological properties is not available), but for shallow sources one should be aware of the issues outlined above. In such cases the inversion should be carefully performed with more constraints (constrained inversions if there is any clue about the nature of the source, or possibly the inclusion of tilt) or in much lower frequency band, less sensitive to structural heterogeneities.
(7) Lowest misfit: The lowest misfit values are not synonym for the best solution. The lowest residuals obtained for the shallow source for model S1 shows the highest discrepancy between the original and retrieved STFs. This is an important remark and caution should be exercised when interpreting results. The summary of our synthetic tests outlined above suggests that the unconstrained inversions for shallow sources with approximate velocity models cannot guarantee correct source solutions. However, it is important to mention that this strong conclusion is obtained from testing inversions in a standard LP frequency band (0.2-2 Hz). It is likely that much better results could be obtained by shifting the scope toward much longer periods and including into inversions both translational ground motion and tilt. We investigated the sensitivity of the moment-tensor inversion solution to the choice of velocity models under a volcanic context. Four models with increasing geological complexity have been used in our synthetic test. Since both the source location and inversion are jointly affected by the uncertainties in the velocity model, we suggest that, when possible, LP events should be first located by other location methods (such as amplitude decay, cross-correlation coefficient and semblance, for example Cannata et al. 2013 and references therein), and then inverted for the best source position. In this way, computational resources for calculating many Green's functions could be attributed for testing various velocity models. If it is not possible to carry out the extensive synthetic testing and there are clues about the nature of the source mechanism, we suggest performing a constrained MT inversion in order to find out the most plausible source mechanism. This is in agreement with the suggestions given in Lokmer et al. (2007) and Bean et al. (2008). Solutions obtained for a shallow source and a homogenous model (S1) tend to overestimate the real amplitude of the STF (in both constrained and unconstrained inversions), hence caution should be exercised when estimating the gas/fluid volumes (possibly) involved in the generation of LP events. Under the presence of shallow unconsolidated volcanic materials, especially when considering shallow sources, all the tested velocity models delivered high uncertainties in the results. In particular, the solution including forces (MT + F) led to an incorrect source mechanism, so we suggest to avoid this type of inversions. Although the results obtained by MT-only inversion exhibit significantly less fluctuation and smaller departure from the true mechanism, they often include a large amount of spurious shear component. Based on the outlined observations, we propose to use the MT-only inversion for deep sources, while the constrained inversion (possibly combined with MT-only) should be used for shallow sources.
It is important to remind that our tests are performed in a typical LP frequency band (0.2-2 Hz), so our conclusions are applicable to this frequency band only. A recent research presented the case where there exists energy in the very-low frequency band (f < 0.1 Hz) of some LP events (Thun et al. 2015). In such a cases, the strategy may be extending the band of inversion towards the low frequencies, which are less sensitive to the structural heterogeneities. Also, apart from the translational signals, the tilt could be included into inversions (e.g. van Driel et al. 2015 andMaeda et al. 2011). Such a joint translational/rotational very-low frequency inversions appear to be a necessary step towards improving our ability to reliably determine a source mechanism from recorded data, and it will be the subject of future work.
Further works on the understanding of the material properties of volcanoes and their response to waves with wider frequency content would strongly improve our understanding of the physical mechanism beyond LP events generation.