*4.3. Acoustic Multipole Sources and Infrasound Waveform Inversion*

The increasing availability of infrasound data over the last decade has allowed important advances in volcano acoustics, and helped overcoming some of the limitations of past studies based on single-station infrasound records. While true amplitude attenuation cannot fully be accounted for without resorting to numerical methods [44,45], using signals from multiple sensors across a network can help mitigating the bias introduced by station location. If available, multiple signals can be simultaneously inverted to constrain volume or mass flow rate source time functions. Johnson et al. [46] performed multi-station inversion solving separately for pure monopole, and pure 3D dipole sources using data from a three-station deployment at Mt. Erebus. They suggested that a more realistic acoustic radiation would be modelled as the combination of the sound produced by the combination a volumetric source (monopole) and a directional force field (dipole); however, inversion for such a source would require data from at least four stations, which were not available. More recently, infrasound waveform inversion was performed considering multipole sources. A second-order multipole source (i.e., combined monopole and dipole) is well-justified when short-duration acoustic signals associated with non-turbulent plumes are analysed. Higher-order source terms (i.e., quadrupoles) should be considered when infrasound is recorded during sustained eruptions, which may generate continuous jet noise due to internal turbulence within the fluid flow. Studies of the sound radiated by sustained volcanic eruptions, until now, have focused on the characterisation of their spectral signatures in analogy to similarity spectra of aeroacoustic jet-noise, or on investigation of the asymmetrical distribution of the recorded amplitudes [1]. We are not aware of any studies, to date, that performed waveform inversions to constrain third-order acoustic multipole sources, thus including a quadrupole term.

A pressure time series, *p*, can be written as the linear superposition of three contributions (Equation (7)), i.e., a monopole, a horizontal dipole and a vertical dipole. The vertical dipole term is often neglected in infrasound studies; vertical dipoles cannot be reliably constrained by data gathered exclusively on the ground, and their inclusion would introduce errors in source estimates [34]. Kim et al. [34] and De Angelis et al. [47] neglected the vertical dipole component in infrasound waveform inversions at Tunghuraua volcano (Ecuador) and Santiaguito (Guatemala); they combined Equation (7) with Equations (4) and (5) to obtain:

$$p(\mathbf{r},t) = \frac{1}{2\pi r} \left[ \dot{S} \left( t - \frac{r}{c} \right) + \frac{x}{cr} \mathcal{D}\_x \left( t - \frac{r}{c} \right) + \frac{y}{cr} \mathcal{D}\_y \left( t - \frac{r}{c} \right) \right] \tag{10}$$

which represents the acoustic pressure radiated in a half-space (a homogeneous and isotropic atmosphere with constant sound speed, bounded by a flat topography) by a combined monopolehorizontal dipole, and recorded at distance *r* from the source. A system of linear equations can be obtained by discretising Equation (10):

$$p\_i^k = \frac{1}{2\pi r\_i} \left[ m\_1^k + \frac{\varkappa\_i}{cr\_i} m\_2^k + \frac{y\_i}{cr\_i} m\_3^k \right] \tag{11}$$

where *p<sup>k</sup> <sup>i</sup>* is the *k*th value of pressure recorded at the *i*th station with coordinates (*xi*, *yi*), at distance *ri* from the source; *m<sup>k</sup> <sup>i</sup>* = [*m<sup>k</sup>* <sup>1</sup>, *<sup>m</sup><sup>k</sup>* <sup>2</sup>, *<sup>m</sup><sup>k</sup>* <sup>3</sup>]=[*S*˙, *<sup>D</sup>*¨ *<sup>x</sup>*, *<sup>D</sup>*¨ *<sup>y</sup>*] is a vector of parameters that represent monopole and dipole strengths. Equation (11) can be written in matrix form:

$$\mathbf{P}^k = \mathbf{G}\mathbf{m}^k\tag{12}$$

where **P***<sup>k</sup>* is a vector of *n* values of pressure measured at *n* stations. Equation (12) can be iteratively solved over all samples in each of the *n* time-series, for the model vector **m***k*:

$$\mathbf{m}^k = \mathbf{G}^{-1} \mathbf{P}^k \tag{13}$$

Equations (10)–(13) provide a framework for inversion of acoustic infrasound waveforms for a multipole source. In Figure 7, we show an example of multipole source inversion for an explosion recorded at Fuego (Guatemala) in May 2018. The time history of monopole (i.e., atmospheric volume flow rate) and dipole (i.e., the magnitude of the dipole force vector) strengths are shown in Figure 7a; it should be noted that the vector *m<sup>k</sup>* <sup>1</sup> retrieved by iteratively solving Equation (13), provides a time history of the rate of atmospheric mass displacement, which is then translated into a volume flow rate using the atmospheric density at vent elevation. This formulation assumes that the volume of displaced atmosphere is equivalent to the volume of material ejected from the vent. In Figure 7b, we show the dominant dipole direction over the source time function. Figure 7c shows the fit between the inversion model and observations. The quality of fit is measured, following De Angelis et al. [47], using variance reduction:

$$VR = \frac{1}{N} \sum\_{i=1}^{N} \left[ 1 - \frac{\int (s\_i - s\_i)^2}{\int (d\_i)^2} \right] \tag{14}$$

where *si* and *di* are the synthetic and data at each of the *N* stations, and integrals are performed over the duration of the signals. For comparison, in Figure 7d,e, we show the results for a monopole-only source, obtained by excluding any dipole terms in Equation (11) from the inversion. Comparison of volume flow rates obtained from multipole (black) and monopole-only (yellow) inversions (Figure 7d,e) shows that the inclusion of dipole terms does not produce significant differences in volume flow rate, although it improves waveform fit between model and data.

**Figure 7.** Comparison of acoustic multipole and monopole-only source inversions for an explosion recorded in May 2018 at Fuego, Guatemala: (**a**) monopole and (horizontal) dipole source time functions obtained from multipole inversion; (**b**) dominant direction of the horizontal dipole over the source time function (grey shaded segment in (**a**)); (**c**) waveform fit between data and model; (**d**) comparison of volume flow rates obtained from monopole-only (yellow) and multipole (black) inversions; and (**e**) waveform fit between data and model for monopole-only inversion. Note that these inversions are based on a representation of the acoustic wavefield as in Equation (10). VR is variance reduction (please, see Equation (14) in main text).

All methods discussed thus far to evaluate eruption source parameters from infrasound data do not consider the effects of topography and local atmospheric conditions on the propagation of the acoustic wavefield. Waveform inversion as described in Equations (10)–(13) assumes half-space sound propagation, and relies on the simplest possible 1D formulation of the atmosphere's GF, that is one that only depends on the distance between source and receiver. Kim and Lees [48] and Lacanna and Ripepe [44] first demonstrated by means of numerical modelling the effects of both crater morphology and local volcano topography on infrasound radiation. More recently, Kim et al. [45] introduced an acoustic waveform inversion method that exploits a numerical approximation to the 3D atmospheric GFs computed by a Finite Difference Time Domain (FTDT) scheme. A GF represents the impulse response of the atmosphere, which is the acoustic pressure generated by an impulsive source, and recorded at a given distance from it. GFs are calculated by solving a set of first-order, velocity-pressure coupled differential equations [49], accounting for the propagation of the acoustic wavefield over realistic topography. In these models at local scale (less than about 10 km source-receiver distance), the atmosphere is considered homogeneous and non-moving, with constant sound speed [45]. Figure 8 illustrates results of FDTD modelling of infrasound propagation at Fuego (Guatemala); both the sound pressure level (i.e., acoustic wave intensity measured in decibels; SPL in Figure 8a), and the GFs computed at different locations on the volcano (Figure 8b) demonstrate complex scattering and attenuation effects. Numerical modelling, thus, reveals that infrasound waveform complexity may, at least partially, reflect topography rather than source effects. Kim et al. [45] and Fee et al. [50] included 3D GFs into a new acoustic source inversion workflow, akin to ordinary seismic moment tensor inversion in seismology. A time series of acoustic pressures, *p*, produced by a monopole and recorded at given location from the source can be written as a convolutional model:

$$p(\mathbf{r}, t) = G(\mathbf{r}, t; \mathbf{r}\_0, t) \* S(t) \tag{15}$$

where *S*(*t*) is the source mass flow rate (i.e., monopole strength) and *G*(**r**, *t*;**r0**, *t*) is the Green's function excited by a source at position **r0** and evaluated at position **r**. Equation (15) can be written in matrix form:

$$\mathbf{d} = \mathbf{Gm} \tag{16}$$

where **d** is a vector of recorded pressures at different stations (locations), **G** is a matrix of GFs at each of these locations, and **m** is a vector of the unknown monopole source strength. Equation (16) can be solved for **m** using an iterative least square solver [51]. In Figure 9a,b, we show an application of monopole inversion using 3D, numerical, GFs at Fuego volcano (Guatemala), and compare the results with those from traditional monopole inversion using 1D GFs. The main observation is that volume flow rates calculated without considering topography are overestimated by about 130% when compared to inversion using 3D GFs. Similar results were reported by Kim et al. [45] at Sakurajima volcano, Japan, where overestimates of volume flow rate were up to 155%.

Numerical simulations of acoustic infrasound over realistic topography are rapidly becoming commonplace due to rapid advances in computational methods, and increasing availability of High Performance Computing resources at comparatively low cost. Several ongoing studies are investigating multipole source inversion using 3D GFs (e.g., [52,53]), and successful attempts to include a vertical dipole component have been made [52]. Including dipole terms in the inversion appears to marginally reduce estimates of fluid flow rates, and has variable influence on the quality of waveform fit between data and the final inversion model. The influence of station density and azimuthal coverage is a factor likely to play an important role on the stability of multipole source inversions although the still relatively small number of experiments with high enough quality datasets means that debate on this subject remains open. For example, due its intrinsic directivity, dipole radiation can only be detected at certain locations, making network azimuthal coverage vital to correctly resolve its orientation and strength. The extent to which the contribution of dipole terms can be confidently resolved during inversion, also depends on the quality and resolution of the digital topography used in GFs numerical modelling, in particular in the near-vent region where more pronounced terrain discontinuities and gradients are encountered. Frequent cloud cover and rapid changes in crater morphology make the availability of up-to-date, high-resolution, digital terrain models difficult, posing challenges to the implementation of acoustic source inversion in (near) real time. Preliminary results from both Iezzi et al. [52] and Diaz-Moreno et al. [53], however, suggest that for impulsive and short-duration explosions, if the effects of topography are correctly accounted for, monopole-only inversions are stable and may be a reasonable approximation to the infrasound source mechanism. Finally, until now, numerical models of acoustic wave propagation have only considered 1D homogeneous representations of the atmosphere with no wind or density changes. This is a reasonably well-justified assumption in the very near-field, although its validity may break down as source-receiver distances approach several kilometres.

**Figure 8.** Numerical simulation of acoustic wavefield propagation at Fuego volcano (Guatemala): (**a**) distribution of Sound Pressure Level (i.e., acoustic wavefield intensity measured in decibels relative to background atmospheric pressure) around the source (white star); and (**b**) 3D Green's Functions excited by a monopole located at the source position (white star) and recorded at different stations across a local infrasound network (station positions on the volcano indicated in (**a**)).

**Figure 9.** Monopole source inversion using 3D numerical Green's functions for an explosion recorded at Fuego volcano (Guatemala): (**a**) comparison of volume flow rates retrieved by monopole source inversion using 3D and 1D Green's functions (black and yellow, respectively); and (**b**) plot of waveform fit between recorded infrasound and best model from 3D monopole inversion. VR is variance reduction (please, see Equation (14) in main text).

We recommend that future infrasound studies should be supported by data from high-density networks, and exploit better digital terrain models to: (i) resolve the relative influence of monopole and dipole terms in evaluating volcanic emissions; and (ii) assess the importance of considering the

3D nature of dipoles on the results and stability of inversions. We suggest that the influence of variable atmospheric conditions should be further explored, as well. The theory and methods of acoustic inversion need to be extended beyond case studies of short-duration transients, to include the more complex infrasound associated with sustained eruptions for which internal plume turbulence may have significant influence on the radiated acoustic wavefield. Finally, we speculate that, despite several caveats, real-time implementation of acoustic source inversion in an operational sense is within close reach. Similar to seismic moment tensor inversion, libraries of GFs could be implemented and stored for use in (near) real-time inversions, and periodically updated if required. We stress the temporal resolution offered by acoustic infrasound in retrieving eruption source parameters in order to inform ash plume rise and transport models remains unmatched, and can play a crucial role for rapid assessment of airborne eruption hazards.

#### **5. Conclusions**

We have presented a comprehensive overview of the most recent developments in the use of acoustic infrasound to assess volcanic emissions, and discussed the theoretical framework of linear acoustic wave theory and its application to volcanic explosions. A wealth of studies over the past two decades have used acoustic infrasound data ranging from single-station pressure time-series to high-density, multi-station datasets to decipher the complexity of the processes that radiate sound at volcanoes. Fluid flow velocities and volume and mass flow rates can be retrieved by means of inversion based on equivalent representations of the acoustic source mechanisms as multipole series. We have discussed the recent introduction of numerical modelling to approximate the atmosphere's impulse response in the presence of realistic topography, and how it has been integrated into inversion workflows. A large body of literature demonstrates the potential of acoustic infrasound to provide rapid estimates of eruption source parameters and their use as input into models of ash plume rise. The still unmatched temporal resolution offered by infrasound makes its use in (near) real time attractive for rapid assessment of airborne hazards during volcanic crises.

**Author Contributions:** All authors contributed equally to the ideas and data analyses presented in the manuscript. S.D.A. initially wrote the manuscript; and all authors worked on final revisions and copy-editing.

**Funding:** Silvio De Angelis and Alejandro Diaz-Moreno are funded by NERC grant number NE/P00105X/1. Luciano Zuccarello has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 798480.

**Acknowledgments:** The authors thank D. Fee and A. Iezzi (Geophysical Institute at the University of Alaska, USA), M. Ripepe and G. Lacanna (Dipartimento di Scienze della Terra, University of Firenze), and J.B. Johnson (Boise State University) for many insightful discussions during the preparation of this manuscript. This manuscript was partly written during a period of research leave spent by S. De Angelis at the University of Alaska Fairbanks; S. De Angelis is indebted to UAF for their generous hospitality and support during his stay. All data used in this manuscript are publicly available through the facilities of the IRIS Data Management Center (http://ds.iris.edu/ds/nodes/ dmc/data/#requests) with the exception of those in: (1) Figure 2b,c,e,f available via direct request to the Istituto Nazionale di Geofisica e Vulcanologia, INGV, Sezione di Catania; and (2) Figure 1b, available via direct request to the authors, or publicly from the National Geoscience Data Centre, UK (https://www.bgs.ac.uk/services/ngdc/). The IRIS DMC is funded through the Seismological Facilities for the Advancement of Geoscience and EarthScope (SAGE), Proposal of the National Science Foundation under Cooperative Agreement EAR-1261681.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
