A nonlinear approach to assess network performance for moment-tensor studies of long-period signals in volcanic settings

SUMMARY

The characterization of seismic source mechanisms and geometries is of critical importance for understanding the underlying causes and physical processes of low-frequency seismicity at volcanoes. Because observational data are often limited in volcanic environments by logistical constraints, we use synthetic modelling to investigate the capability of seismic networks to properly resolve source mechanisms. For 16 synthetic networks with as many as 40 stations, and variable distance and azimuthal distributions, we perform nonlinear moment-tensor inversion for six input source models. Using a grid search for source type and constrained waveform inversions provides a quantitative measure of source mechanism reliability. If the source location is assumed to be correct, results suggest that complete azimuthal coverage, with stations located at different distance ranges from the source will allow for a higher resolution recovery of the source-time function. In general, a similar degree of uncertainty characterizes configurations with as many as 40 stations and as few as eight stations. Although the level of uncertainty in the source-time function increases when fewer than eight stations are used in the inversions, sources are still recoverable when as few as four stations are used. Deviations from this general trend are present across the different input source models, with performance of configurations with fewer than eight stations showing a strong dependence on the source type.

Inverse theory, Waveform inversion, Volcano seismology

1 INTRODUCTION

Volcanic eruptions involve a complex interaction of multiphase fluids with solid conduit walls. Seismology has shown great promise for modelling conduit dynamics (e.g. Chouet 1996), as seismic signals can carry information related to fluid transport phenomena, such as coalescence of bubbles, ascent of gas slugs or sealing and pressurization phenomena. Source modelling of seismic events beneath volcanoes can therefore produce models that correspond to fluid pathway processes and geometries, leading to a better understanding of the underlying causes and physical processes connected to eruption activity. A clearer picture of the dynamics of magmatic systems is achieved by the analysis of low-frequency events, which include tremor, long-period (LP, 0.2 – 2 s) and very-long-period (VLP, 2 to 100 s) events. Full-waveform moment-tensor inversion of seismic data has been a principal technique used to successfully derive seismic source mechanisms at a number of active volcanoes for both VLP and LP events. LP events are typically interpreted as due to pressurized fluids (hydrothermal, gas or magma) in conduits and cracks (e.g. Dawson et al.2011; Chouet & Dawson 2016) although some have argued for a slow, brittle failure mechanism (Eyre et al.2015). We refer to Chouet & Matoza (2013) for a recent review.

However, moment-tensor inversions are non-unique and bear intrinsic uncertainties due to simplifying assumptions about the source and structure such as inaccuracies in the velocity models and the dependence on network configuration. In order to advance source modelling, it is important to quantify the sensitivity of model results to these controlling variables. Various studies have focused on the influence of near-surface volcanic structure on moment-tensor inversions and demonstrated that overly simplified velocity structure can have a detrimental effect on source inversions (e.g. Bean et al.2008; Trovato et al.2016), especially when the source is shallow (Trovato et al.2016). When the effects of near-surface structures are not taken into account, they can lead to the emergence of incorrectly oriented source geometries, spurious single forces and incorrect source-time functions. Incorrect solutions may also be related to volcano topography. Topographical scattering can produce path effects which can leak into the source solutions (O'Brien & Bean 2009). A priori structural information and topography are therefore essential for modelling the nature of shallow short-duration, volcano-seismic signals, such as LPs. In addition to near-surface low-velocity structure and topography, network configuration has a great influence on the solution of LP and moment-tensor inversions. The capability of seismic networks to accurately resolve source mechanisms and locations is one of the primary factors we need to consider in order to develop an accurate interpretation of seismic sources. The effect of volcano seismic network configurations on LP source mechanisms has not been thoroughly investigated to date but is especially important to understand in volcanic environments where accessibility and equipment are often limited.

We use a nonlinear inversion for seismic source types that allows for comparison of, for example, double-couple (DC) or compensated linear vector dipole (CLVD) sources versus isotropic components and provides a quantitative description of moment-tensor uncertainty (Waite & Lanza 2016). The importance of showing the variations of the misfits over the space of source types is also highlighted in Alvizuri & Tape (2016) where uncertainty estimation of moment tensors of small magnitude events at Uturuncu volcano, Bolivia, is represented by the variance reduction for the misfit function between the observed and synthetic waveforms.

In this study, we focus on the sensitivity of source inversion results to network configuration using a variety of idealized networks, a known velocity model, and the topography of Pacaya volcano, Guatemala. Pacaya volcano is the subject of a parallel study in which the nonlinear procedure described in this paper is applied to the source inversion of explosion-related long-period events (see Lanza & Waite 2018, in press). We use synthetic modelling to explore the network configuration effects on moment-tensor solutions. We test six different synthetic input source models with as many as 40 stations at variable azimuths to assess the influence of the network parameters (i.e. number of stations and azimuthal coverage) that can contribute to an increase/decrease on the accuracy of the source-type solution. For each network configuration we perform constrained inversions and search over all possible moment-tensor types and orientations at the source model input location and compare the model fits for different station configurations.

2 INVERSION METHODOLOGY

2.1 Full-waveform inversion procedure

We followed an inversion approach that is similar to other studies (see e.g. Auger et al.2006; Waite et al.2008; Dawson et al.2011; Richardson & Waite 2013) where the inversion is performed in the frequency domain to reduce the computation time and permits a grid search over a large volume. The representation theorem to describe the displacement field generated by a point source in the frequency domain is

\begin{eqnarray*} {u_n}\ \left( {\vec{r},\omega } \right) = {F_{p\ }}\ \left( \omega \right) \cdot {G_{np}}\left( {\vec{r},\omega } \right) + \ {M_{pq}}\left( \omega \right) \cdot {G_{np,q}}\left( {\vec{r},\omega } \right), \end{eqnarray*}

(1)

where |${u_n}( {\vec{r},\ \omega } )$| is the Fourier transform of the n component of the seismogram, p and q are direction indices x, y and z, M_pq(⁠|$\omega $|⁠) is the Fourier transform of the time history of the pq-component of the moment tensor and |${G_{np}}( {\vec{r},\ \omega } )$| is the Green's function that relates the n-component of displacement at the receiver position, |$\vec{r}$|⁠, with the moment at the source position and |$\omega $| is angular frequency. The notation ,q indicates spatial differentiation with respect to the q-coordinate.

In matrix form, eq. (1) becomes

\begin{eqnarray*} {\boldsymbol{U\ }}\left( \omega \right) = {\boldsymbol{G}}\left( \omega \right){\boldsymbol{s}}\left( \omega \right), \end{eqnarray*}

(2)

where U is the N_t x 1 vector of Fourier-transformed ground displacement components, G is the N_t x 9 (or 6) matrix of Fourier transforms of the Green's functions, s is the 9 (or 6) x 1 vector of Fourier transformed force and moment-tensor components (or moment-tensor components only) and N_t is the number of observed seismic traces.

We solved for the full spectra of the bandpass filtered (0.5 – 4 Hz) seismograms and invert at each frequency separately by minimizing the least-squares problem:

\begin{eqnarray*} {\boldsymbol{s\ }}\left( \omega \right) = {\left[ {{\boldsymbol{G}}{{\left( \omega \right)}^T}{\boldsymbol{G}}\left( \omega \right)} \right]^{ - 1}}\ {{\boldsymbol{G}}^T}\left( \omega \right){\boldsymbol{U}}\left( \omega \right). \end{eqnarray*}

(3)

We then applied the inverse Fourier transforms to obtain the time estimates of both the source components and synthetic seismograms. In this paper, we refer to ‘unconstrained’ inversion when we invert for different combinations of the nine model parameters (six moment-tensor components and three single forces), using eq. (3) without applying any a priori constraints to the solution. The term ‘constrained’ inversion is used to describe inversion in which the moment tensor for the source-time function is fixed. Constrained inversions are performed as part of the nonlinear inversion procedure which involves a grid search over all possible moment-tensor types and orientations. Details of the constrained inversion are given in the next section.

For both unconstrained and constrained inversions, synthetic Green's functions were computed with the 3-D finite-difference method of Ohminato & Chouet (1997). We used a cosine smoothing function to derive the Green's functions and to assure stability in the inversion. This wavelet has a time-constant period of 1 s (1 Hz) to approximate the peak frequency of the LP events recorded at Pacaya volcano. The Green's functions convolved with the cosine function represent the elementary source-time functions used in the inversion. We used a model that includes the 3-D topography of Pacaya, derived from a digital elevation map from 2006 with a resolution of 10 m. The model domain is centred on the active summit crater of the volcano and has lateral dimensions of 4 km by 4 km and a vertical extent of 3.5 km. This yields a model with 401 × 401 × 351 nodes spaced 10 m apart. All station locations were rounded to the nearest node, and topography was resampled to match the node grid. The node spacing satisfies the criterion of minimum number of grids per wavelength of 25 established by Ohminato & Chouet (1997). The model is also wide enough to minimize edge reflections of the boundaries while including all the stations.

We use a homogeneous velocity model with P-wave velocity of 978 m s⁻¹, S-wave velocity of 565 m s⁻¹ and density of 1750 kg m⁻³. These values are based on the mean of the first five layers (first 130 m in depth) of the shear wave velocity structure derived from Love and Rayleigh dispersion curves (Lanza et al.2016). We acknowledge that a homogeneous model is not correct, and, as demonstrated by Trovato et al. (2016) at Etna volcano, while a simple velocity model could be used for deep sources, for shallow sources the best available tomographic model should be considered, and more constraints should be applied. However, because this study is concerned primarily with source model sensitivity to station configuration, we deem a simple velocity model appropriate.

To test the ability to find the input location for each of the point sources, we conducted a grid search over a volume of 740 m x 740 m x 500 m centred on the summit vent. In order to reduce the total number of points inverted for, we skipped every other node so that the effective spacing in the search volume was 20 m. Even with the 20 m grid spacing, and excluding points above the topography, the number of possible point source locations becomes quite large, reaching a total of 22 088 sources. If for each source location we calculate nine Green's functions, corresponding to the nine independent model parameters (six moment-tensor components and three single forces), we are required to generate 198 792 Green's functions. Following Chouet et al. (2005), we used the reciprocal relation (Aki & Richards 2002) between source and receiver, which greatly reduces the number of finite-difference runs required to derive the Green's functions. Treating each station as a point source and each potential source as a receiver reduced the sources to 40 (equal to the highest number of stations used in the inversions) and left us with 22 088 receivers. For each source we computed reciprocal Green's functions for each of three orthogonal single forces, a total of only 120 runs with the finite-difference code. In order to compute the full set of moment-tensor Green's functions we calculated the spatial derivative between appropriate single force components in adjacent nodes in the potential source volume. We validated the Green's functions calculated this way by running forward models for selected nodes at the corners (bottom and top) of our point-source volume and comparing the results to those obtained by the reciprocal method. We found that the results from the two methods were indistinguishable for the eight nodes that we tested.

2.2 Nonlinear inversion procedure

The unconstrained moment-tensor inversion at hundreds of points within a 3-D model provides a spatial estimate of the location and its corresponding uncertainty. Moment-tensor orientation and mechanism type can be obtained at the best-fitting source location using point-by-point eigenvector decomposition of the source-time function (e.g. Chouet et al.2003). The moment tensor can then be interpreted in terms of one or more physical components of the source model, such as cracks, pipes or spherical sources. While this method has been very successful, it does not provide a means to assess the range of acceptable source types. We seek to expand on this approach by exploring the uncertainty in the source type through a grid search over all possible moment-tensor types and orientations (see also Waite & Lanza 2016). In this nonlinear inversion, the single forces are kept as free parameters.

The search over source types is determined using the fundamental lune source-type definition of Tape & Tape (2012). Similar to a Hudson-type plot (Hudson et al.1989), the lune maps the distribution of all possible moment-tensor source types onto a surface. Source types are defined by two parameters, the spherical coordinates γ and δ, which are defined by the ratios of the moment-tensor eigenvalues. The longitude parameter, γ, ranges from –30° to 30°, and it is defined as

\begin{eqnarray*} {\gamma} = \mathrm{ arctan} \left( {\frac{{{{ - \ }}{{{L}}_{{1}}}{{ + \ 2}}{{{L}}_{{{2\ }}}}{{ - \ }}{{{L}}_{{3}}}}}{{\sqrt {{3}} {\boldsymbol{\ }}\left( {{{{L}}_{{1}}}{\boldsymbol{\ }}{{ - \ }}{{{L}}_{{3}}}} \right)}}} \right), \end{eqnarray*}

(4)

where |${L_1}$|⁠, |${L_2}$| and |${L_3}$| are the eigenvalues of each point of the source-time function.

The latitude parameter, δ, ranges from –90° to 90°, and it is derived from the co-latitude parameter, β, for which

\begin{eqnarray*} {\beta} =\mathrm{ arccos }\left( {\frac{{{{{L}}_{{1}}}{{ + \ }}{{{L}}_{{2}}}{{ + \ }}{{{L}}_{{3}}}}}{{\sqrt {{3}} {{\ |}}\left| {{{\boldsymbol \Lambda}}} \right|{{|}}}}} \right){{\ \ \ \ \ {\rm with}\ \ \ }}\left| {\left| {{{\boldsymbol \Lambda}}} \right|} \right|{{ = }}\sqrt {{{{L}}_{{1}}}^{{2}}{{ + }}{{{L}}_{{2}}}^{{2}}{{ + }}{{{L}}_{{3}}}^{{2}}} \end{eqnarray*}

(5)

\begin{eqnarray*} {\rm \delta} = 90^{\circ} {\rm - \beta }, \end{eqnarray*}

(6)

where |${L_1}$|⁠, |${L_2}$| and |${L_3}$| are the eigenvalues of each point of the source-time function. The latitude δ is zero for deviatoric patterns, and both δ and γ are zero for DC mechanisms.

The search over γ and δ uses the surface spline method described in Tape & Tape (2012), which effectively reduces the number of points required to evenly sample the distribution over an evenly spaced grid. We evaluate the full range of γ, but because the lower half of the lune is simply the opposite sign of the upper half (e.g. positive isotropic decrease versus negative), we evaluate δ from 0° to 90°.

In addition to exploring the full source-type space, the range of possible orientations for each source is also examined. To do this, we rotate the moment tensor at 10° intervals using a sequence of three rotations about the initial coordinate system of the moment tensor. Complete sampling of the symmetric tensor orientations requires exploring three rotations: a full 360° range about one axis, for example, the z-axis, a range from 0° to 180° of dip of z to z′, and finally, 90° of rotation about the new, rotated z′ axis (Goldstein et al.2002; Waite & Lanza 2016). The 10° interval was chosen for computational reasons and after having assessed, through multiple tests with finer intervals, that although there are small changes in the misfit values there were no significant changes in the pattern of misfits. This involves 5832 combinations of rotation angles combined with 223 γ–δ pairs, for a total of 1 300 536 trial moment-tensor solutions.

For each individual moment source type in the gridded lune we invert for the source-time function that provides the best fit to the data by fixing the moment tensor. We refer to Waite & Lanza (2016) for the details of the inversion approach.

2.3 Evaluation of the results

Quantification of uncertainties in seismic source inversions are challenging. Although recently full probabilistic approaches via Bayesian inference have emerged (Stähler & Sigloch 2014), it is still not well understood how to best quantify error for moment-tensor solutions. Here we follow a common methodology adopted in many applications of moment-tensor inversion at volcanoes: the overall selection of the best solution is based on a weighted squared error between the observed (in our case synthetic) and reconstructed data, that is either normalized by station or channel (Ohminato et al.1998; Chouet et al.2003). In this study we adopted the squared error measure described as E₁ where the large amplitude signals on the near-source stations dominate the error so that mismatches between low amplitude data and synthetics do not dominate the error:

\begin{eqnarray*} {{\rm{E}}_{\rm{1}}}{\rm{ = }}\frac{{\mathop \sum \nolimits_{{\rm{\mathit{ n} = 1}}}^{{{\rm{N}}_{\rm{\mathit{ r}}}}} \mathop \sum \nolimits_{{\rm{\mathit{ p} = 1}}}^{{{\rm{N}}_{\rm{\mathit{ s}}}}} {{\left( {{\rm{u}}_{\rm{\mathrm{\mathit{ n}}\mathit{ }}}^{\rm{0}}\left( {{\rm{\mathit{ p}}}\Delta {\rm{t}}} \right){\rm{ - u}}_{\rm{\mathit{ n}}}^{\rm{s}}\left( {{\rm{\mathit{ p}}}\Delta {\rm{t}}} \right)} \right)}^{\rm{2}}}}}{{\mathop \sum \nolimits_{{\rm{\mathit{ n} = 1}}}^{{{\rm{N}}_{\rm{\mathit{ r}}}}} \mathop \sum \nolimits_{{\rm{\mathit{ p} = 1}}}^{{{\rm{N}}_{\rm{\mathit{ s}}}}} {{\left( {{\rm{u}}_{\rm{\mathit{ n}}}^{\rm{0}}\left( {{\rm{\mathit{ p}}}\Delta {\rm{t}}} \right)} \right)}^{\rm{2}}}}}, \end{eqnarray*}

(7)

where |$u_n^0(p\Delta t)$| is the pth sample of the nth data trace, |$u_n^s(p\Delta t)$| is the pth sample of the nth synthetic trace, N_s is the number of samples in each trace and N_r is the number of three-component receivers.

For the unconstrained inversion results, the influence of the number of free model parameters in each source model was evaluated by calculating the Akaike's information criterion (AIC; Akaike 1974), expressed as

\begin{eqnarray*} \mathrm{ AIC} = {N_\mathrm{ t}}\ {N_\mathrm{ s}}\mathrm{ ln}{E_\mathrm{ 1}} + 2{N_\mathrm{ m}}{N_\mathrm{ f}}, \end{eqnarray*}

(8)

where E₁ is the squared error from eq. (7), |${N_\mathrm{ m}}$| is the number of source mechanisms and |${N_\mathrm{ f}}$| is the number of frequencies in the passband of interest. In general, the use of more free parameters is considered justified when both the squared error and AIC are reduced.

3 SYNTHETIC MODELLING

Small networks of temporary or permanent stations are common on active volcanoes, but it is not clear what the ideal network geometry or the recommended number of stations is for this type of study. Thus, we conducted an evaluation of the ability of ‘synthetic’ networks with as many as 40 stations to resolve key components of the moment tensor. This approach is similar to the work by Dawson et al. (2011) at Augustine volcano, Alaska, but in our case, we perform a more exhaustive exploration of possible moment-tensor source types. Importantly, this study does not consider unknowns in the velocity model which can have important effects on waveforms in the LP band, especially at distances greater than a few km (e.g. Richardson & Waite 2013) and for sources at shallow depths (Trovato et al.2016). For longer period (VLP) sources, the velocity errors are insignificant given the much longer wavelengths.

We computed a simulated wavefield for six input source models: (1) a crack dipping 10° with dipole magnitude ratios of [1, 1, 2]; (2) a vertical crack striking 80° E with dipole magnitude ratios of [1, 3, 1]; (3) an isotropic source [1, 1, 1]; (4) a pure CLVD [2, −1, −1]; (5) a linear vector dipole (LVD) [1, 1, 0] and (6) a DC [1, 0, −1] (in an east-north-up coordinate system). For each model, Green's functions were computed with and without the addition of small single force components in x, y and z. The amplitude of the single forces was two orders of magnitude smaller than the moment components.

To synthesize a realistic set of seismograms, the Green's functions calculated for each input source model were convolved with a portion of the M_xx component of the source-time function obtained from the unconstrained inversion of data collected at Pacaya (see Lanza & Waite 2018, in press). Synthetics were computed for an input location situated 260 m directly beneath the summit vent of Pacaya volcano, at 2330 m a.s.l. The choice of the depth of the input location was functional to have all of the considered source points below topography, while still considering the case for a fairly shallow source. To simulate noisy traces in real data, white noise bandpass filtered in the LP band was added to the final synthetic waveforms. The noise amplitude was calculated using the average of the root mean squared of the synthetic signal amplitudes from the four closest stations. In this way, the noise amplitude was constant across the network.

We explored the effect of network configuration for sixteen subsets of the synthetic network, each consisting of three to 40 stations with variable azimuthal distribution, and whose geometries are shown in Fig. 1. Configuration C8, indicated by red triangles in Fig. 1, was chosen to be identical to the actual local network deployed at Pacaya volcano in 2013 (see Lanza & Waite 2018, in press). Although we inverted for 16 different subnets for each source input model for a total of 96 configurations, we show here the detailed analysis and results for six representative models for clarity. A summary of the results of the unconstrained inversion from all the configurations can be found in the Supporting Information of this paper (Table S1 and Fig. S1). In the Supporting Information we also include an example of waveform fits for the unconstrained moment tensor only and moment tensor with single forces inversions for configuration C8 and the six input models (Fig. S2, Supporting Information) and the results obtained from the constrained inversions for all configurations and the six input models (Figs S3–S8, Supporting Information).

Figure 1.

Subnet geometries of the 16 station configurations considered for the synthetic modelling. Input source model location is represented by the black star (central inset). Red triangles indicate the position of C8 which is identical to the network deployed at Pacaya during the 2013 October–November field campaign.

Although all configurations were useful for exploring different combinations of number of stations and azimuthal coverage as summarized in Fig. 6; for readability and clarity, in this paper we focus on subnets labelled in Fig. 1 as C1, C4, C5, C6, C8 and C10 as they are representative of complete or poor azimuthal coverage and span a variable number of possible stations employed. Subnet C1 and C4 have both the same number of stations and the same azimuthal coverage, but differ in the distance from which the stations are located with respect to the source; subnet C5 has 16 stations at varying distance from the source with a complete azimuthal coverage; subnet C6 has 16 stations clustered north and south the input source model; subnet C8 has only four stations with fairly good azimuth coverage, and subnet C10 includes all 40 stations. Subnet C8 reflects the location of the temporary seismic network deployed at Pacaya.

3.1 Unconstrained inversion results

All 96 configurations are inverted both with and without single forces in the solution, here referred as MT-only inversion, and MT + F inversion. For each of the selected subnets of stations, the residual errors and locations of the best-fitting unconstrained inversions (E_1min solutions) are reported in Table 1. Table 1 shows that the best-fitting locations of all trials for both MT-only and MT+F inversions cluster to a small area around the original input model location, within approximately 20 m of the synthetic source location in both north–south and east–west directions. A comprehensive table and a summary figure with the results from all 96 configurations are available in the Supporting Information of this paper (Table S1 and Fig. S1). Larger variations, up to 180 m, are observed for depth and generally for those configurations that have fewer stations or a lower azimuthal coverage as shown by configurations C13–C16 in Table S1 of the Supporting Information. For these configurations higher location misfits are predominantly associated with the MT+F inversions.

Table 1.

Synthetic inversion results for the unconstrained inversions.

			MT-only solutions (MT-only synthetics)							MT +F solutions (MT+F synthetics)
Subnets	# of stations	Input source	E1	AIC	E–W position relative to input source [m]	N–S position relative to input source [m]	Depth position relative to input source [m]†	Eigenvalue ratios* (unconstrained inversion)	Eigenvalue ratios* (constrained inversion)	E1	AIC	E–W position relative to input source [m]	N–S position relative to input source [m]	Depth position relative to input source [m]†	Eigenvalue ratios* (unconstrained inversion)	Eigenvalue ratios* (constrained inversion)
1	10	Isotropic (1, 1, 1)	0.069	−14 810	0	0	20	[1 0.98 0.95]	[1 1 1]	0.046	−15 897	0	20 N	−40	[1 0.98 0.88]	[1 1 1]
		Crack (1, 1, 2)	0.050	−16 976	0	0	0	[2 1.01 1]	[2 1.14 1.1]	0.042	−16 608	0	0	0	[2 0.9 0.89]	[2 1.14 1.1]
		Crack (1, 3, 1)	0.042	−18 116	0	0	0	[3 1 0.95]	[3 0.84 0.76]	0.035	−17 824	0	0	0	[3 0.96 0.74]	[3 0.84 0.76]
		LVD (1, 1, 0)	0.040	−18 420	0	0	20	[1 0.99 0.15]	[1 0.99 −0.16]	0.034	−18 039	0	0	0	[1 0.98 0.53]	[1 0.97 −0.73]
		CLVD (2, −1, −1)	0.040	−18 333	0	0	0	[0.46 −0.45 −2]	[1.16 1.12 −2]	0.035	−17 773	0	0	0	[0.20 0.17 −2]	[1.16 1.12 −2]
		DC (1, 0, −1)	0.033	−19 691	0	0	0	[0.68 0.15 −1]	[1 0 −1]	0.029	−19 079	0	0	0	[0.33 0.27 −1]	[1 0 −1]
4	10	Isotropic (1, 1, 1)	0.289	−5220	0	0	40	[1 0.97 0.94]	[1 1 1]	0.186	−6603	0	0	20	[1 0.78 0.72]	[1 1 1]
		Crack (1, 1, 2)	0.165	−8964	0	0	20	[2 0.97 1.03]	[2 1.14 1.1]	0.142	−8420	0	0	20	[2 1.03 0.99]	[2 1.14 1.1]
		Crack (1, 3, 1)	0.296	−5058	0	20 N	0	[3 1.28 1.1]	[3 0.84 0.76]	0.252	−4577	0	20 N	0	[3 1.39 0.94]	[3 0.84 0.76]
		LVD (1, 1, 0)	0.219	−7067	0	0	20	[1 0.95 0.09]	[1 0.96 0.03]	0.176	−6978	0	0	20	[1 0.86 0.53]	[1 0.96 0.03]
		CLVD (2, −1, −1)	0.215	−7177	0	0	20	[1.12 0.96 −2]	[1.09 0.91 −2]	0.175	−7012	0	20 N	20	[0.65 0.57 −2]	[1.09 0.91 −2]
		DC (1, 0, −1)	0.204	−7533	0	0	0	[0.70 0.20 −1]	[1 0 −1]	0.177	−6930	0	0	−20	[0.54 0.38 −1]	[1 0 −1]
5	16	Isotropic (1, 1, 1)	0.119	−19 662	0	0	0	[1 0.98 0.95]	[1 1 1]	0.085	−21 764	0	20 N	−20	[1 0.86 0.81]	[1 1 1]
		Crack (1, 1, 2)	0.084	−23 405	0	0	−20	[2 1.03 1]	[2 1.14 1.1]	0.077	−22 723	0	0	−20	[2 1 0.98]	[2 1.14 1.1]
		Crack (1, 3, 1)	0.097	−21 788	0	0	0	[3 0.98 0.88]	[3 0.84 0.76]	0.088	−21 322	0	0	0	[3 0.89 0.79]	[3 0.84 0.76]
		LVD (1, 1, 0)	0.097	−21 791	0	0	0	[1 0.98 0.03]	[1 0.99 −0.16]	0.081	−22 243	0	0	0	[1 0.87 0.36]	[1 0.96 0.03]
		CLVD (2, −1, −1)	0.096	−21 936	0	0	0	[1.01 0.91 −2]	[1.16 1.12 −2]	0.087	−21 433	0	0	0	[0.38 0.37 −2]	[1.16 1.12 −2]
		DC (1, 0, −1)	0.080	−23 869	0	0	0	[0.83 0.09 −1]	[1 0 −1]	0.075	−22 988	0	0	0	[0.37 0.24 −1]	[1 0 −1]
6	16	Isotropic (1, 1, 1)	0.119	−19 638	0	20 N	0	[1 0.99 0.95]	[1 0.92 0.84]	0.079	−22 488	0	20 N	−20	[1 0.86 0.82]	[1 1 1]
		Crack (1, 1, 2)	0.086	−23 096	0	0	−20	[2 1.03 1]	[2 1.14 1.1]	0.078	−22 660	0	0	−20	[2 0.99 0.97]	[2 1.14 1.1]
		Crack (1, 3, 1)	0.073	−24 925	0	0	0	[3 1 0.84]	[3 0.84 0.76]	0.065	−24 551	0	0	0	[3 1.15 1.07]	[3 0.84 0.76]
		LVD (1, 1, 0)	0.099	−21 617	0	0	0	[1 0.98 0.03]	[1 0.99 −0.16]	0.084	−21 865	0	0	0	[1 0.87 0.33]	[1 0.99 −0.16]
		CLVD (2, −1, −1)	0.094	−22 158	0	0	0	[1.02 0.92 −2]	[1.09 0.91 −2]	0.086	−21 554	0	0	0	[0.33 0.33 −2]	[1.09 0.91 −2]
		DC (1, 0, −1)	0.053	−28 231	0	0	0	[0.76 0.12 −1]	[1 0 −1]	0.049	−27 552	0	0	0	[0.39 0.24 −1]	[1 0 −1]
8	4	Isotropic (1, 1, 1)	0.024	−6921	0	0	−20	[1 0.98 0.95]	[1 1 1]	0.011	−7467	0	20 N	−80	[1 0.96 0.92]	[1 1 1]
		Crack (1, 1, 2)	0.017	−7723	0	0	0	[2 1.06 1.01]	[2 1.14 1.1]	0.009	−7831	0	20 S	20	[2 1.13 1.09]	[2 1.14 1.1]
		Crack (1, 3, 1)	0.011	−8978	0	0	0	[3 0.94 0.87]	[3 0.84 0.76]	0.005	−9331	0	0	0	[3 0.88 0.52]	[3 0.84 0.76]
		LVD (1, 1, 0)	0.020	−7387	0	0	−20	[1 0.95 0.15]	[1 0.99 −0.16]	0.007	−8460	0	0	40	[1 0.80 0.65]	[1 0.99 −0.16]
		CLVD (2, −1, −1)	0.019	−7487	0	0	0	[1.11 1.01 −2]	[1.16 1.12 −2]	0.008	−8236	0	0	20	[0.21 0.13 −2]	[1.16 1.12 −2]
		DC (1, 0, −1)	0.010	−9241	0	0	0	[0.81 0.05 −1]	[1 0 −1]	0.005	−9399	0	0	0	[0.45 0.10 −1]	[1 0 −1]
10	40	Isotropic (1, 1, 1)	0.134	−50 640	0	0	0	[1 0.98 0.95]	[1 1 1]	0.100	−56 890	0	0	−20	[1 0.87 0.83]	[1 1 1]
		Crack (1, 1, 2)	0.091	−60 914	0	0	0	[2 1.02 1]	[2 1.14 1.1]	0.086	−60 780	0	0	−20	[2 0.98 0.96]	[2 1.14 1.1]
		Crack (1, 3, 1)	0.118	−53 938	0	0	0	[3 0.99 0.94]	[3 0.84 0.76]	0.113	−53 647	0	0	0	[3 1.08 1]	[3 0.84 0.76]
		LVD (1, 1, 0)	0.110	−55 740	0	0	0	[1 0.99 0.05]	[1 0.99 −0.16]	0.099	−57190	0	0	0	[1.00 0.87 0.24]	[1 0.99 −0.16]
		CLVD (2, −1, −1)	0.105	−57 028	0	0	0	[1.03 0.92 −2]	[1.16 1.12 −2]	0.100	−56 888	0	0	0	[0.47 0.44 −2]	[1.16 1.12 −2]
		DC (1, 0, −1)	0.092	−60 651	0	0	0	[0.81 0.09 −1]	[1 0 −1]	0.089	−59 885	0	0	0	[0.42 0.24 −1]	[1 0 −1]

^*The eigenvalues ratios are always associated with the solution at the input source location.

^†Positive depth values indicate shallower depth with respect to the input location.

For comparison purposes and given the similarities in the best-fitting centroid location to the input source location for all source models and subsets, we evaluated the consistency of the source-time function at the input source location instead of the best-fitting centroid obtained from the unconstrained inversion. The consistency of the source-time function over its duration is measured through point-by-point eigenvector decomposition (Chouet et al.2003; Waite et al.2008). The ratios of the eigenvalues of the moment tensor give insight into the source mechanism while the eigenvectors give the orientation of the principal axes. We only consider points with amplitudes of at least 80 per cent of the peak amplitude for this analysis. For the selected configurations (C1, C4, C5, C6, C8 and C10) the median ratios of minimum and intermediate eigenvalues to the largest are reported in Table 1. Note that the eigenvalue ratios refer to the solutions at the input source location rather than at each best-fitting centroid (E_1min solutions), and that we also report the eigenvalue ratios derived from the constrained solutions (see Section 3.2 for details on the constrained inversion). For the MT-only solutions, the ratios of the eigenvalues are in agreement with the input source mechanisms. On the other hand, for some configurations and source models (e.g. C4 for DC and CLVD source models) higher deviations from the eigenvalues of the input models characterize the MT+F solutions, even though the misfits are lower. In order to statistically determine the utility of extra model parameters, we calculated AIC which is reported in Table 1 (and Table S1 of the Supporting Information for all configurations). We observe that both the squared error and the AIC values are generally reduced for isotropic sources and for all source models when the number of stations is low (i.e. configuration C8). This behaviour reflects how the addition of more free parameters when using a limited number of stations is more likely to provide a numerically stable solution by absorbing model errors (Bean et al.2008; De Barros et al.2013). The apparent better fits are however contrasted by the recovery of a source mechanism that clearly has a worse fit with respect to the input model as shown by the higher deviations from the eigenvalues. This leads to an important conclusion: the smallest inversion misfit or lowest AIC does not always correspond to the most accurate source mechanism (as in Matoza et al.2015; Trovato et al.2016).

We further examined the stability of source-time history by calculating the lune longitude γ and latitude δ for each point of the source-time history of each moment tensor. The use of γ–δ representation can be very useful to quantify the stability of the derived moment-tensor source type. We represent each point of the source-time history in the γ–δ space where points with negative δ were projected to the upper portion of the plot (Fig. 2). In general, MT+F models produce a lower misfit, but a worse match to the input source model as indicated by the eigenvalue ratios (Table 1) and the γ–δ distribution on the lune (Fig. 2).

Figure 2.

(a–f) Point-by-point decomposition of the source-time functions represented in the γ–δ plane for the results of the unconstrained inversion for each input source model. For clarity, only the six most representative station configurations are shown [subnets C1 (a), C4 (b), C5(c), C6 (e), C8 (d), and C10 (f)]. Colours refer to the six input source types considered in the study; while solid circles represent MT-only solutions, triangles indicate MT+F solutions. Points with negative δ are projected to the upper portion of the plot. Source-time functions for both MT-only and MT+F solutions are also shown for each configuration and input source model. Black lines indicate MT-only solutions, while grey lines represent MT+F solutions.

One reason behind the better recovery for MT-only inversions, even when the synthetic model contains nine components, could be attributable to the small magnitude of the single forces compared to the moment-tensor components. In order to better understand the role of the single forces, we performed additional inversions with synthetics that included single force components of higher magnitude. These tests indicate that the source mechanisms are very difficult to retrieve if the single forces are in phase with the moment-tensor components.

Although it has been demonstrated that the single forces may absorb model errors while providing a numerically stable solution (Bean et al.2008; De Barros et al.2013), in reality, as shown here and by recent findings from Trovato et al. (2016), single forces can increase the discrepancy of the retrieved from the true solution. In fact, Trovato et al. (2016) suggest to generally avoid the MT + F inversions. While LP sources in nature may indeed involve single forces, we find that they may be very difficult to model. However, we caution that our treatment of the single forces was somewhat simplified and perhaps unrealistic. For example, we would not expect single forces for DC or CLVD mechanisms. But since our principal concern in this study is the influence on network configuration on model resolution, we chose to keep the synthetic single forces the same across all the models for consistency.

The point-by-point eigenvalue decomposition of the MT-only moment-tensor source-time functions and the γ–δ representation shows that all source models are well recovered with the best fits for the isotropic source and poorer fits for the CLVD and DC sources. Greater inconsistency, indicated by greater scattering in both γ and δ values, is associated with subnets that either have an uneven azimuthal coverage (C6), or fewer stations (C8).

If the azimuthal coverage and the number of stations is the same, increasing the source–receiver distance generally has a detrimental effect on the stability of the moment-tensor solution (e.g. subnets C1 through C4 in Figs S3–S8 of the Supporting Information). Velocity structures at volcanoes can commonly be complicated by near-surface low-velocity layers that introduce errors in velocity models and require consideration of signal contamination from path effects and with real data a higher degradation effect should be expected. Other studies have therefore suggested deploying stations as close to the source as possible (e.g. Neuberg & Pointer 2000). Decreasing the number of stations from 10 to 8 (subnets C11 and C12), while maintaining good azimuth coverage does not decrease the resolution of the source-time functions. This pattern is common for all six source input models.

One of the advantages of synthetic modelling is that there are no uncertainties in the velocity model, meaning the differences in the source-time functions across the different subnets must be attributed to other factors. Figs 2(a)–(f) show the time history of the source on the moment components for the six selected configurations for all six input sources, where the only parameter that is changed in the inversions is the network configuration. In the next section, we further investigated these differences in the source mechanisms when considering different subnets by performing nonlinear inversions for source type.

3.2 Constrained inversion results

We next perform a grid search using the constrained inversion approach described above. The unconstrained inversions over a volume of points showed that there is very little uncertainty in the best-fitting source locations. For this reason, and for the ease of comparison with the unconstrained inversion results, we used the input model location to investigate the constraint on the source type. We perform the constrained inversions using only the moment components because of the difficulty in resolving single forces noted in Section 3.1. The full grid search over all possible source types was conducted for each of the subnets and source models.

In Fig. 3 we show the error for each point on the lune source type plot for the selected subnets; for each moment-tensor type the orientation that gives the lowest misfit is used. The point-by-point source-time function eigenvector analysis for the MT-only unconstrained inversions that were plotted in Fig. 2 are shown in the appropriate plots with white circles for comparison. The contours define the misfit values for all moment-tensor orientations for each γ–δ pair in the lune plot. The lack of stability in the unconstrained inversion for subnet C6 is reflected in the lack of a resolved moment tensor in the constrained inversion for almost all input source models.

Figure 3.

Misfit by fixed moment-tensor solution, plotted together with the point-by-point mechanism type for the unconstrained inversion as white circles. Only the six subnets discussed in the text are shown here for clarity. Note that for the CLVD input source the white circles, representing the point-by-point decomposition for the unconstrained inversion, are plotted on the left-hand side of the lune as a result of the projection of the negative points in the upper-half portion of the lune.

To give a sense of the waveform fits on individual channels, we computed cross-correlation coefficients and lag times between the synthetic input data, and the modelled waveforms from both the unconstrained inversion, and the constrained inversion at the lowest misfit value of the contour error maps shown in Fig. 3. Fig. 4 shows an example of the waveform fits for network configurations C5, C6, C8 and C10 with the CLVD input source. For clarity, we show only station number 14 as this station is common to all subnets considered, except for subnet C8, in which a different set of stations are used as it reflects the actual local network installed at Pacaya volcano. However, station 14 is located only few meters away from the location of station P03 of the seismic array, therefore we deemed acceptable to use station P03 for subnet C8 (see inset in Fig. S2 of the Supporting Information for station location).

Figure 4.

Data (black line), unconstrained solutions (red line) and best-fitting constrained solutions (blue line), are shown for a CLVD source for each of the four configurations discussed in the paper: C5, C6, C8 and C10. The unconstrained inversion results are in good agreement with the best-fitting constrained inversion models.

The cross-correlation coefficients for both the unconstrained and constrained inversion are above 0.95 with no lag time for the higher amplitude vertical and north components (Table 2). The east component exhibits a significantly lower correlation and lag times ranging from −0.2 to 0.4 s. Station P03 is located to the north of the source, and for this mechanism the east component has little energy. This indicates a high level of agreement between the waveform fits from the unconstrained and best fixed inversions.

Table 2.

Correlations between data from one common station (S14) and synthetic waveforms from the unconstrained and constrained inversion models, and corresponding lag times, for each channel used.

Inversion type	Subnet	Channel	Corr. coef.	Lag [s]	Channel	Corr. coef.	Lag [s]	Channel	Corr. coef.	Lag [s]
unconstrained	C5	S14 Z	0.96	0.00	S14 N	0.96	0.00	S14 E	0.67	0.02
	C6	S14 Z	0.97	0.00	S14 N	0.95	0.00	S14 E	0.61	−0.02
	C8	P03 Z	1.00	0.00	P03 N	1.00	0.00	P03 E	0.97	0.00
	C10	S14 Z	0.96	0.00	S14 N	0.95	0.00	S14 E	0.62	0.02
constrained (min E1min)	C5	S14 Z	0.96	0.00	S14 N	0.95	0.00	S14 E	0.58	0.02
	C6	S14 Z	0.96	0.00	S14 N	0.95	0.00	S14 E	0.59	0.00
	C8	P03 Z	0.99	0.00	P03 N	1	0.00	P03 E	0.89	0.02
	C10	S14 Z	0.96	0.00	S14 N	0.95	0.00	S14 E	0.58	0.04

The nonlinear analysis confirms the general observations from the unconstrained inversions. Clearly, subnet C6, in which we have an uneven station distribution with respect to azimuthal coverage, shows fairly broad uncertainties (Fig. 3). When considering networks with similar azimuthal distribution but different number of stations, similar fits are obtained for subnet C5 (16 stations) and subnet C10 (40 stations) indicating that the resolution of the moment tensor is not necessarily improved when a larger number of stations are included in the inversion.

The relative resolving ability of the various networks can be demonstrated quantitatively by examining the area of the lune that falls below a chosen cut-off of E₁. This allows for comparison across all configurations for a given input source and for all sources for a given configuration. To facilitate this, we chose to adopt a cut-off value based on the minimum value of E₁ for each configuration and add a constant value of 0.1 to the minimum. This somewhat arbitrary cut-off value is meant to define a misfit level for which the model no longer fits the data, as the synthetics are too different from the data. However, an appropriate cut-off is subjective as it is also a function of the noise level in the data. For example, a relatively high E₁ for noisy data can produce a fairly good fit, while the same E₁ for very clean data would not match the synthetics as well as in the first case. Fig. 5 shows examples of traces from a high S/N and a low S/N input waveform compared with the corresponding modelled synthetic. Both produce good fits although the noisy trace is associated with a much higher misfit value.

Figure 5.

Cut-off dependence on noise level. Left-hand and right-hand panels show waveform fits comparison between two traces that have different signal-to-noise ratio. Note how the fits are comparable (fairly good fits) in both cases, but the misfit value is much higher in the noisy trace.

The portion (expressed in percentage) of the lune that falls below the error cut-off for each input source and each configuration are shown in Fig. 6(a). Fig. 6(b) is a graphical representation of each point in the crack (1,1,2) panel of Fig. 6(a). Fig. S9 in the Supporting Information shows for the same graphical representation of Fig. 6(b) for all input source models. In general, results computed with larger numbers of stations are better (they have a small lune portion below the chosen cut-off value) than results computed with just a few stations. For example, C8, which has four stations has more dispersed γ–δ pairs and wider range of possible mechanism types as shown in Fig. 3 with lune portions of about double the average of the other configurations with more stations (Fig. 6a). Higher uncertainties and poorly constrained source-time functions always characterize configurations with fewer than four stations.

Figure 6.

(a) Portion of the lune space that falls below the E₁cut-off for all 16 configurations across the six input source models for the constrained solutions. The size of the circles is proportional to the number of stations of each network and colours denote maximum azimuthal gap for the configuration. (b) Graphical representation of the portion of the lune space below E₁ cut-off for the crack (1,1,2) input source model (second panel from the left in a).

Fig. 6(a) also highlights the importance of azimuthal distribution. As expected, configurations with good azimuthal coverage (e.g. C5 or C7 which have maximum azimuthal gaps between adjacent stations <130°), perform much better than configurations with the same number of stations but larger azimuthal gaps (e.g. configuration C6). For all mechanisms, a similar percentage is calculated for configurations that include as many as 40 stations and as low as eight stations. This is good news for those cases when deployments are limited due to logistical issues or lack of resources. However, for configurations within this range in terms of number of stations, but that have a higher percentage of stations located far from the source, the misfit values are higher and show higher uncertainties (e.g. Fig. 6a for configuration C4), while the γ–δ solution is fairly stable and well resolved (Fig. S5, Supporting Information). This, once again, underlines the subjectivity of the cut-off value and its dependence on noise level in the data, and how not only the misfit value but also the misfit pattern along with the γ–δ pairs is to be considered to assess the overall quality of the solution.

From Fig. 6(a) we also observe how the quality of the model fits is strongly dependent on the input source model. For example, a crack mechanism, across all configurations, is recovered and resolved with a better accuracy (lower percentage below the cut-off, Fig. 6a) than source types with a non-isotropic component (e.g. a DC or CLVD model). For those source models that do not have an isotropic component, the source inversion generally has a wider range of acceptable models than for those with an isotropic component (Fig. 6a; Figs S7 and S8, Supporting Information). The nonlinear inversions for CLVD and DC input models suggest a wider range of possible mechanism types fit the data across all configurations, even though the waveform fits at the lowest misfit values show high correlation (Table 2). For example, the range of acceptable solutions for a CLVD input source model includes sources in the crack/isotropic regions in some cases. As might be expected, this problem is most significant for network configurations with smaller numbers of stations and poor azimuthal coverage, but it is present in the other configurations as well.

Finally, we explored the effect of source depth on the resolution of moment-tensor solutions for the different synthetic station configurations. This might be especially important in the case of volcanic sources where all stations are on the edifice at a range of elevations. We applied the same inversion procedure to (1) deeper input source models located at 1400 m a.s.l., about 250 m below the lowest station and (2) shallower input source models located at 2330 a.s.l. The shallowest point is 260 m below the highest point in the topography model. Previous work has shown that shallow sources can be difficult to resolve, especially when near-surface low-velocity layers are present (Trovato et al.2016). In our tests, the results show similar findings to the conclusions of Trovato et al. (2016). The deeper source shows lower uncertainties when the misfit is evaluated in a source-type space (Fig. 7). Because topography and velocity models are known, this behaviour at shallower depths may be attributed to an incomplete coverage of the focal sphere that ultimately leads to greater uncertainty in the mechanism. For example, C4 has 10 stations, all at about the elevation of the deeper source, but they are all well below the shallow source. On the other hand, C1 has 10 stations all at about the elevation of the shallow source. For the CLVD source (Fig. 7b), this is reflected in the opposite behaviour of E₁ for C1 and C4. The crack mechanism does not show the same reversal in C1 and C4 and that is partly due to the fact that a crack-oriented E–W is not so dependent on the distribution of points on the focal sphere.

Figure 7.

Percentage of the lune misfit below the chosen cut-off for a crack and CLVD input sources for the two different tested depths at 2330 (grey) and 1400 m (red).

The nonlinear inversion procedure not only allows us to quantify the uncertainty of the source type, but also provides a means to explore the range of acceptable source orientations. We analysed the relation between misfits and orientation of the source mechanisms for the best-fitting γ–δ pair for all station configurations and for each of the source input models, with the exclusion of the isotropic case. Perhaps unsurprisingly we found that the subnets with a better moment-tensor resolution appear to have a better control on the geometry of the source, but overall, we found that the orientations of the moment tensors are very well constrained. Fig. 8 shows misfit plots for the moment-tensor orientations in the case of a crack input model (1,3,1) together with the rose diagrams of φ (azimuth) and ϑ (inclination angle) for moment-tensor solutions whose orientations have E₁ values below the misfit cut-off chosen for each network configuration. We show the orientation of the maximum, intermediate and minimum eigenvectors to the right of the rose diagrams. We observe that the variations in misfit for moment orientations for C6 and C8, where the stations are either poorly distributed or in small number, are significantly broader with a wider range of solutions that present low misfits. Moreover, the farther the stations (C4), the higher are the misfit values, although, the crack orientation is still well recovered.

Figure 8.

(a) Misfit plot of the crack model (1,3,1) orientations for the six selected network configurations. (b) Rose diagrams for azimuth (φ) and inclination angle (ϑ) of moment-tensor solution with E₁ below the chosen cut-off value. Maximum, intermediate and minimum eigenvectors are shown as blue, green and yellow bars, respectively.

In general, this suggests that errors not only in the source type, but also in the geometry of the retrieved mechanism, are more likely introduced when azimuth and/or number of stations is inadequate.

4 CONCLUDING REMARKS

We investigated the capability and accuracy of limited seismic networks that are common on volcanoes to recover LP source mechanisms using synthetic modelling. We examined six mechanism types with 16 network configurations using a nonlinear inversion approach. Because we explore the entire solution space, this inversion approach for source type provides quantitative constraints on source model uncertainty, which can aid in the interpretation of mechanism types.

One observation from this study is that larger numbers of stations are nearly always better, but when network sizes are equal, those distant stations relative to the source can degrade model resolution (e.g. subnet C1 versus C4, Fig. 2). In our synthetic tests the degradation of the source recovery is minor, but we expect to have a greater influence on the moment-tensor solution in real cases where errors in the velocity model due to the complicated near-surface velocity structure can lead to source contamination from path effects that increase with distance from the source.

Another observation regards the network azimuthal distribution. When azimuthal coverage is poor, that is, when the gap is greater than 130°, the model uncertainty increases dramatically. Moreover, solutions were generally good across all input source models when more than eight stations where considered, but no significant improvements were reached by employing more than 16 stations.

Uncertainty, as measured by the percentage of the total model space below a chosen misfit value increases gradually with decreasing numbers of stations, but then usually rises when just three stations are used. Based on these observations, networks with as few as four stations should be capable of modelling moment tensors provided they are well distributed around the source. However, we find that the resolution of the recovered source-time function is dependent on the source type; lower uncertainties are generally associated with source mechanisms having a larger isotropic component; CLVD sources are particularly difficult to constrain. Nonetheless, the waveform fits for the best-fitting constrained inversion model agree well with the unconstrained inversion solution.

Our results also show that the lowest misfit values are not necessarily representative of the best solutions, as in Bean et al. (2008) and Trovato et al. (2016). The inclusion of single forces may increase the discrepancy between source model and true solution, especially if they are in phase with the moment components, a finding also in agreement with suggestions by Trovato et al. (2016). However, contrary to the findings of spurious DC components in MT-only inversions noted by Trovato et al. (2016), we observe a predominant contamination of an isotropic component for DC and CLVD solutions in which a wider range of solutions can even approach the crack domain on the lune.

In this work, the use of known velocities and densities in a model with realistic 3-D topography allows for isolation of model sensitivity to network geometry. Previous work has examined the important role of near-surface low-velocity layers in resolving sources in the LP band. For example, Trovato et al. (2016) explored velocity models of varied complexity and observed that major uncertainties occur for shallower sources even with a high number (16) of stations. They conclude that unconstrained inversions with approximate velocity models cannot guarantee a correct solution and suggest constrained inversions should be used for shallow sources. Clearly, the errors associated with unmodelled structure will be magnified in cases with smaller numbers of stations, and especially networks with more distant stations. Given that there will always be unmodelled structure, we suggest that those variables that can be precisely controlled, such as network configuration, and station placement, should be explored and optimized when possible. This study demonstrates that constrained inversions and the mapping of the fit of the seismic moment tensor in source-type space serves as a guide to planning station location on volcanoes and later interpretation of the results.

SUPPORTING INFORMATION

Table S1. Synthetic inversion results for the unconstrained inversions.

Figure S1. E₁ misfit values for MT-only and MT+F unconstrained inversions of all 96 tested configurations. Blue and magenta circles indicate that the minimum solution occurred at the input source location for MT-only and MT+F inversion, respectively.

Figure S2. Waveform fits for configuration C8 between the six input models (black) and the synthetics computed as a result of the unconstrained inversion at the input source location (red).

Figure S3. Misfit by fixed moment-tensor solution for the isotropic input model (1, 1, 1), plotted together with the point-by-point mechanism type for the free inversion as white circles. Labels in the upper right-hand corner of each lune indicate the station configuration number. We refer to Fig. 1 in the paper for subnet geometries.

Figure S4. Misfit by fixed moment-tensor solution for the horizontal crack input model (1, 1, 2), plotted together with the point-by-point mechanism type for the free inversion as white circles. Labels in the upper right-hand corner of each lune indicate the station configuration number. We refer to Fig. 1 in the paper for subnet geometries.

Figure S5. Misfit by fixed moment-tensor solution for the vertical crack input model (1, 3, 1), plotted together with the point-by-point mechanism type for the free inversion as white circles. Labels in the upper right-hand corner of each lune indicate the station configuration number. We refer to Fig. 1 in the paper for subnet geometries.

Figure S6. Misfit by fixed moment-tensor solution for the LVD input model (1, 1, 0), plotted together with the point-by-point mechanism type for the free inversion as white circles. Labels in the upper right-hand corner of each lune indicate the station configuration number. We refer to Fig. 1 in the paper for subnet geometries.

Figure S7. Misfit by fixed moment-tensor solution for the CLVD input model (2, −1, −1), plotted together with the point-by-point mechanism type for the free inversion as white circles. Labels in the upper right-hand corner of each lune indicate the station configuration number. We refer to Fig. 1 in the paper for subnet geometries.

Figure S8. Misfit by fixed moment-tensor solution for the DC input model (1, 0, −1), plotted together with the point-by-point mechanism type for the free inversion as white circles. Labels in the upper right-hand corner of each lune indicate the station configuration number. We refer to Fig. 1 in the paper for subnet geometries.

Figure S9. Graphical representation of the portion of the lune space below E₁ cut-off for all input source models and network configurations.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

This work was supported by National Science Foundation Award No. 1053794. The manuscript has greatly benefitted from comments by Philip Dawson and two anonymous reviewers.