Recent Developments and Applications of Acoustic Infrasound to Monitor Volcanic Emissions

: Volcanic ash is a well-known hazard to population, infrastructure, and commercial and civil aviation. Early assessment of the parameters that control the development and evolution of volcanic plumes is crucial to effective risk mitigation. Acoustic infrasound is a ground-based remote sensing technique—increasingly popular in the past two decades—that allows rapid estimates of eruption source parameters, including ﬂuid ﬂow velocities and volume ﬂow rates of erupted material. The rate at which material is ejected from volcanic vents during eruptions, is one of the main inputs into models of atmospheric ash transport used to dispatch aviation warnings during eruptive crises. During explosive activity at volcanoes, the injection of hot gas-laden pyroclasts into the atmosphere generates acoustic waves that are recorded at local, regional and global scale. Within the framework of linear acoustic theory, infrasound sources can be modelled as multipole series, and acoustic pressure waveforms can be inverted to obtain the time history of volume ﬂow at the vent. Here, we review near-ﬁeld (<10 km from the vent) linear acoustic wave theory and its applications to the assessment of eruption source parameters. We evaluate recent advances in volcano infrasound modelling and inversion, and comment on the advantages and current limitations of these methods. We review published case studies from different volcanoes and show applications to new data that provide a benchmark for future acoustic infrasound studies.


Introduction
Approximately 10% of the Earth's population live under the direct threat of one of 1508 active volcanoes [1,2]. Understanding the nature and impact of volcanic hazards and developing monitoring systems to detect and track eruptions in (near) real time, is central to effective early warning and successful risk mitigation. The type and severity of hazards posed by volcanic activity depend on the style of eruption. Effusive eruptions often affect settlements and infrastructure such as during the recent 2018 rift eruption of Kilauea (Hawaii, USA) [3]. At the other end of the spectrum of volcanism, explosive and sustained activity can result in loss of life, such as during the 2018 eruption of Volcan de Fuego (Guatemala), and have a severe impact on local communities [4]. In general, long-lived eruptions can have significant, negative, repercussions on the economy at the national and international scale. During volcanic explosions fragments of pulverised rock, referred to as volcanic ash, are ejected from one or more active vents at high speed, driven into the atmosphere by hot gasses. It is well-known that airborne volcanic ash presents a direct threat to aviation; it can cause disruption to flight operations such as during the 2010 eruption of Eyafjallajokull (Iceland) and damage to infrastructure, with consequent loss of revenue [5].
During volcanic crises nine international Volcanic Ash Advisory Centers (VAAC), part of the International Civil Aviation Organization's (ICAO) Air Navigation Commission, are deputed to dispatch Volcanic Ash Advisories (i.e., alerts on the presence of airborne volcanic ash and estimates of its movement) to aviation authorities [6]. Numerical models are routinely used by VAACs to forecast atmospheric dispersal of ash clouds in order to evaluate the potential impacts of eruptions on aviation. The volume of erupted material, the rate at which mass is ejected from volcanic vents and the initial height of eruption plumes, collectively referred to as eruption source parameters [5], are key inputs into these models. Assessment of these parameters relies on empirical laws, which were derived from measurements of well-documented eruptions [7][8][9][10]. These empirical relationships are integrated in routine protocols to estimate mass eruption rates (MER) from measurements of ash plume height [6], despite large uncertainties when results are compared with values obtained from direct measurements of eruption deposits and durations [5]. Mastin et al. [5] discussed whether MER could be alternatively estimated by means of numerical modelling: based on the results of 1D numerical simulations, they argued that empirical equations still perform better than numerical models in reproducing field observations.
In recent times, acoustic infrasound has emerged as an increasingly popular tool for volcano remote sensing, and its potential to assess volcanic emissions in (near) real time has been extensively investigated. The term infrasound identifies atmospheric acoustic waves with frequencies typically <20 Hz, below the audible range of humans. Volcanoes are prolific radiators of infrasound, generated by large-scale eruptive processes, which inject gas and pyroclasts into the atmosphere causing its rapid acceleration [11]; these low-frequency acoustic waves can travel distances of up to several thousands of kilometres [1] lending themselves to volcano monitoring at different scales, from local to global. The use of acoustic infrasound in regional and local eruption monitoring, with applications in volcano early warning has become increasingly popular [12][13][14][15][16]. A substantial body of literature demonstrate the use of infrasound to detect, locate and track explosive volcanic eruptions, and its potential to provide estimates of eruption source parameters to inform volcanic plume rise and ash dispersal modelling [17][18][19][20][21][22][23].
In this manuscript, we provide a review of past work in the field of volcano infrasound, and describe recent developments of its use to assess eruption source parameters. Here, we focus on applications of acoustic infrasound at local distances (i.e., within 10-15 km from eruptive vents) and their potential to provide robust estimates of the amount of material ejected into the atmosphere during eruptions. We present a summary of the theory of linear acoustics and discuss equivalent models for the representation of volcano-acoustic sources. Methods derived from this theory that exploit volcano-acoustic time series to retrieve eruptive jet velocities, and invert infrasound waveforms to obtain estimates of mass and volume flow rates (of the eruptive mixture of gas and particles) are reviewed. We discuss the recent integration of numerical modelling of acoustic wavefield propagation within these workflows. We show examples from the published literature, as well as applications to new case studies. Finally, we offer a perspective on the potential for use of acoustic infrasound in real-time volcano monitoring.

The Acoustic Fingerprint of Volcanoes
Volcanic activity encompasses a broad spectrum of eruptive styles ranging from effusive discharge of lava onto the Earth's surface, to energetic explosions that inject large amounts of gas-laden pyroclasts in the atmosphere. Acoustic waves are produced when the atmosphere is perturbed from its background state by a moving source causing unsteady fluid motion, and pressure disturbances propagate through the air. Multiple processes occur in volcanic environments making them ideal sources of acoustic waves; these include gas-driven oscillations of the surface of lava lakes or magma columns [24], surface mass movements and impacts (rockfalls, pyroclastic flows and lahars [25]) explosive fragmentation of magma at different scales (from Strombolian to Plinian [17,20]), and vigorous ash-and-gas jetting during sustained eruptions [18]. Volcanoes produce acoustic waves that deliver energy over a broad spectrum of frequencies; however, because of the comparatively large spatial scale of source mechanisms, the majority of sound is radiated in the infrasonic band, approximately between 0.01 and 20 Hz [11]. Infrasound waveforms associated with explosions at volcanoes that exhibit regular activity. Excess pressures are reduced to a distance of 1 km from the vent for comparison across volcanoes: (a) small degassing explosions recorded at Etna, Italy; (b) small degassing explosions recorded at Santiaguito, Guatemala; (c,d) ash-rich impulsive explosions from Tunghurahua, Ecuador, and Fuego, Guatemala; (e,f) large amplitude signals exhibiting blast wave characteristics recorded at Fuego, Guatemala, and Sakurajima, Japan; (g,h) short-duration signals featuring multiple pulses from Fuego, Guatemala; (i) monochromatic signal recorded at Mt. Etna associated with degassing activity from one of the summit vents; abd (j,k) explosion signals with impulsive onsets followed by an extended coda from Fuego, Guatemala, and Sabancaya, Peru, associated with observations of sustained ash plumes.
Volcano infrasound waveforms range from impulsive transients ( Figure 1) to tremor-like sustained signals, with a number of other intermediate types [26]. Waveforms with short duration (5-15 s) and sharp onsets (Figure 1a-f) are frequently generated by impulsive explosion-like sources. These signals may feature compression onsets followed by rarefaction with a nearly symmetrical shape ( Figure 1a); this symmetry is indicative of a flow rate source time function symmetrically distributed, in time, around a peak value [27]. In many instances, waveforms are less symmetrical, exhibiting rapid compression onsets (e.g., Figure 1e) followed by rarefaction phases with reduced amplitudes (Figure 1b-f); waveform asymmetry may represent either a non-symmetrical flow rate source function, or reflect a shock-type source mechanism similar to blast waves produced by chemical explosions [27,28]. Typically, blast waves can be separated from other explosion mechanisms due to their characteristic appearance and much larger peak amplitudes, of the order of several hundreds of Pa within a few hundred metres from the source. Diversity in the characteristics of infrasound signals reflects variety in source mechanisms, and may additionally provide clues on whether explosions are gas-or tephra-rich [29]. Observational evidence suggests that waveforms featuring impulsive onsets followed by several additional pulses (Figure 1g-h), or by a prolonged coda (Figure 1j-k), are frequently associated to tephra-rich plumes with durations up to a few minutes. Open-vent volcanic systems, where active degassing is persistent, often produce nearly monochromatic waveforms [24] due to resonance taking place within the shallow conduit and crater system ( Figure 1i).
Longer duration (several minutes to hours) infrasound is typically associated with sustained open-vent degassing, or continuous eruptive activity. Figure 2 illustrates three examples of extended duration waveforms associated with persistent pulse-like degassing, a sequence of intermittent small Strombolian explosions, and an episode of sustained lava fountaining. A wider variety of tremor-like signals have been reported in the literature; a comprehensive account of these signals is, however, beyond the scope of this manuscript. For additional details, the reader is referred to a number of insightful reviews published in recent years on the nature and characteristics of volcano infrasound [11,26,30,31].

Linear Acoustic Theory and Multipole Acoustic Sources
Fluid flow associated with volcanic explosions causes acceleration of the atmosphere making them efficient radiators of acoustic infrasound. The propagation of density (and thus, pressure) perturbations in a fluid medium at rest can be investigated with the Navier-Stokes equations, by solving an associated wave equation [32]. Early, seminal, work by Lighthill [33] and Woulff and McGetchin [32] set out an elegant formulation for the solution of the acoustic wave equation in free-and half-space, to describe the radiation of sound from monopole, dipole, and quadrupole sources. A monopole is the simplest source of acoustic waves, linked to unsteady fluid flow (i.e., injection of mass) such as during a volcanic explosion, and it is the most efficient form of acoustic radiation. Dipole radiation results from the application of a force field at the source location, and it is obtained as the linear superposition of two adjacent monopoles oscillating out of phase by 180 • . Dipole sources, unlike monopoles, do not introduce net mass at the source; they impart momentum to the fluid causing its movement, and thus, to radiate sound. Woulff and McGetchin [32] suggested that dipole radiation may result from the interaction of fluid flow from volcanic vents with the crater walls. Quadrupole radiation represents turbulence within the gas jet itself, and it results from the combination of two opposite dipoles. Woulff and McGetchin [32], building on previous work by Lighthill [33], and using dimensional analysis, presented a suite of scaling laws that link the acoustic power radiated by monopole (Π m ), dipole (Π d ) and quadrupole (Π q ) sources to gas velocity during fluid flow from volcanic vents: where K m , K d and K q are dimensionless constants, ρ p is the density of the volcanic jet, c is the speed of sound in the atmosphere, v e is the gas exit velocity at the vent, and A v is the cross-sectional area of the vent. Recently, Kim et al. [34] presented a detailed review of the derivation of the formal integral solution to the wave equation in a half-space for monopole and dipole sources, using the Green's Function (GF) solution to the inhomogeneous Helmholtz wave equation [35,36]. Under the assumptions of compact acoustic sources (i.e., their characteristic dimensions are much smaller than the acoustic wavelength) and that far-field volcano infrasound propagates in half-space (i.e., an atmosphere bounded by a flat solid topography), the transient solutions to the wave equation for monopole, and horizontal and vertical dipole radiation at distance r from the source can be written as: where p m , p h , and p v are monopole, horizontal and vertical dipole pressure fields, respectively;Ṡ is the time history of mass acceleration at the source; (x, y, z) are coordinates in a Cartesian system centred at the source location; D x,y,z are dipole momenta (in units of kg m s −1 ) in the (x, y, z) directions; θ is the angle between the axis of the dipole and the vertical direction; z 0 is source elevation; and c is the sound speed. Any infrasound time series can be written as the radiation from a multipole source, that is the linear superposition of a monopole, with horizontal and vertical dipoles (and terms of higher order): The formulation of Equations (4)- (7) represents the foundation for acoustic waveform inversion that is discussed below.

Eruption Source Parameters and Their Potential Use in Ash Plume Modelling
In this section, we discuss the use of infrasound for the evaluation of eruption source parameters, and their integration into models of ash plume rise. We review some case studies and expand on past work to include recent developments in numerical modelling of volcano acoustic wavefields, and their integration into infrasound waveform inversion schemes for the assessment of volume and mass flow rates during explosive eruptions.

Fluid Flow Velocity from Acoustic Infrasound
Early work from Woulff and McGetchin [32] demonstrated the application of Equations (1)-(3) to assess volcanic jet dynamics-in particular fluid flow velocity-from records of the sound emitted by explosions at Stromboli (Italy) and fumaroles at Acatenango (Guatemala). An experiment at Acatenango where acoustic power, the cross-sectional area of the vent, and flow velocity could be directly measured in the field, allowed calibration of the dimensionless constants K m , K d , and K d in Equations (1)-(3), and to infer that the sound radiated in either case had, predominantly, dipolar nature. Several authors [17][18][19][20] built on the original work of Woulff and McGetchin [32] on translating time series of acoustic pressure into gas flow velocities. Acoustic power (Π) can be estimated directly from the recorded infrasound [20] as: where ∆p(t) is a time series of atmospheric pressure, r is the source-receiver distance, ρ atm is atmospheric air density, c is the atmospheric sound speed and τ is the source duration. Once acoustic power is calculated by solving the definite integral in Equation (8), the flow velocity at the source can be easily retrieved by inversion of Equations (1)-(3). Gas velocity is an important eruption parameter easily translated into volume flow rate (for a known cross-sectional area of the vent), which is used as input into one of many empirical models to evaluate the initial height of eruption plumes. The most used of models is the one from Sparks et al. [8]: where Q magma is the magma volume flow rate in units of m 3 /s and H plume is the plume height in units of kilometres. Vergniolle and Caplan-Auerbach [20] used acoustic data recorded during the sub-Plinian phase of the 1999 eruption of Shishaldin volcano (USA) to infer gas velocity and total gas volume erupted; they assumed dipole radiation during the sub-Plinian phase of the eruption, and confirmed that the gas velocities obtained were realistic using independent plume height measurements and Equation (9). They noted that the assumption of dipole radiation would not be appropriate outside the period of sub-Plinian activity; for instance, discrete explosions are more effectively described as monopole sources. A dipole source is often the preferred choice for explosive events lasting up to few tens of seconds associated with plumes rising to their highest elevation within a few minutes. In these situations, the force field exerted by the conduit and crater walls on the volcanic flow, and the interactions of the gas phase with solid pyroclasts are well described by a dipole. During sustained eruptions, lasting up to several hours, infrasound is generated by internal turbulence within the plume, a source best represented as a quadrupole [17].
Caplan-Auerbach et al. [17] also favoured a dipole source to describe the infrasound recorded during the 2006 eruption of Augustine volcano (USA); they analysed nine eruptive events with durations of the order of several minutes that generated largely buoyant plumes reaching elevations of 8-15 km above the vent. Here (Figure 3), we have reproduced their results following the workflow presented in the original manuscript for two of the explosions (Events 1 and 4, in the original publication). Estimates of flow velocity at the vent were obtained, in the original manuscript and, here, by using Equation (8) to calculate acoustic power and then inverting Equation (2). They calculated volume flow rates for all nine explosions assuming a vent with circular cross-section, and used these values to predict plume heights. Modelling the plumes as buoyant thermals rather than continuously fed by eruption (an implicit assumption of Equation (9)), they showed that gas fluxes retrieved from analysis of acoustic data are a good indicator for the rise height of volcanic plumes. Ideally, the choice of a specific source mechanism should be cross-validated using other geophysical measurements, when available. For example, high-resolution visual and thermal infrared imagery can be used to assess eruption dynamics and to calculate the exit velocity of volcanic jets in the near-vent region, two crucial factors in controlling sound generation. Ripepe et al. [18] employed multi-parameter data to investigate the 2011 eruption of Eyjafjallajokull (Iceland); plume heights were retrieved from infrasound-derived fluid flow velocities assuming a dipole acoustic source, and cross-validating these measurements with continuous thermal IR measurements of the plume in the near-vent region. Careful consideration, when using infrasound-derived flow velocities for the assessment of plume rise, should be given to the nature of the plume itself, and whether its initial development is affected by local atmospheric conditions. Empirical relationships such as Equation (9) are usually derived for strong plumes where internal drag is negligible and entrainment of atmospheric air is approximately proportional to the vertical velocity of the flow [37]. For weak plumes, where internal shear is near normal to the plume axis, turbulence may be important and some of the empirical relationships found in the literature may not be appropriate [38]. Furthermore, plumes that develop from low to intermediate velocity volcanic jets are more likely to be affected by local wind conditions in the lower troposphere, causing plumes to bend. These effects are included in more complex numerical models of ash plume rise. Woodhouse et al. [39] developed an integral model of volcanic plume rising in a windy atmosphere and investigated how wind restricted the rise height of volcanic plumes at Eyjafjallajokull. Lamb et al. [19] used the same plume rise model to investigate two events during the 2009 eruption of Mt. Redoubt, USA. Fluid flow velocities were estimated from acoustic infrasound for monopole, dipole and quadrupole sources, and used to produce models of plume rise, including local wind conditions ( Figure 4). Comparison of plume modelling [39] with syn-eruptive Doppler radar measurements of the plumes [40] showed good agreement with the results obtained for infrasound dipole radiation.

Acoustic Monopoles: Volume and Mass Flow Rate
One of the advantages in the use of acoustic infrasound for assessment of volcanic emissions is that its source is directly linked to the unsteady injection of mass into the atmosphere and its acceleration, and to other forces acting in the near-source region. Equations (4)-(6), thus, provide a framework for retrieval of volume or mass flow rates of volcanic material ejected during eruptions. For a pure monopole, the relationship between volume or mass flow rate and the resulting acoustic pressure time series is given by Equation (4), as illustrated in Figure 5; for a source radiating as a monopole in a homogeneous and isotropic atmosphere, Equation (4) can be inverted to retrieve the rate at which atmospheric mass is displaced during an explosion by direct integration over time of the recorded pressure time series. Johnson [41], Johnson and Lees [42], and Johnson and Miller [43] demonstrated the application of the monopole source model to explosions recorded at Erebus (Antarctica), Santiaguito (Guatemala), and Sakurajima (Japan). Direct waveform integration according to Equation (4) rests on the assumption that infrasound amplitudes attenuate proportionally to the inverse of distance from the source (1/r, or geometrical spreading), and neglects the effect of topography on the acoustic wavefield. It is, however, well-known that acoustic wave amplitudes often do not attenuate proportionally to 1/r (e.g., Lacanna and Ripepe [44]), thus limiting the application of the monopole model to single-station recordings. Johnson and Miller [43] at Sakurajima proposed its use with data recorded at sites at proximal distance from the source and with unobstructed line-of-sight to the vent, conditions that provide a reasonable approximation to geometrical spreading and allow discarding major effects of topography. In Figure 6, we show an application of the monopole model to data collected by a temporary infrasound network at Fuego (Guatemala) during a moderate-size explosion in May 2018. This example illustrates well the effects of non-1/r attenuation and topography at different locations around the source. Acoustic waveform shape (Figure 6a-e) varies across the network, evidence of diffraction and scattering from topography; the source time functions of volume flow rate (Figure 6f,j) show variable shapes and inconsistent peak values, suggesting that Equation (4) cannot account for the true attenuation pattern, and for realistic scattering of the acoustic wavefield.

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: 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): where p k i is the kth value of pressure recorded at the ith station with coordinates (x i , y i ), at distance r i from the source; is a vector of parameters that represent monopole and dipole strengths. Equation (11) can be written in matrix form: where P k 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 : 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 k 1 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: where s i and d i 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.  (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: where S(t) is the source mass flow rate (i.e., monopole strength) and G(r, t; r 0 , t) is the Green's function excited by a source at position r 0 and evaluated at position r. Equation (15) can be written in matrix form: 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.  . 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.

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.