Volume Flow Rate Estimation for Small Explosions at Mt. Etna, Italy, From Acoustic Waveform Inversion

Rapid assessment of the volume and the rate at which gas and pyroclasts are injected into the atmosphere during volcanic explosions is key to effective eruption hazard mitigation. Here, we use data from a dense infrasound network deployed in 2017 on Mt. Etna, Italy, to estimate eruptive volume flow rates (VFRs) during small gas‐and‐ash explosions. We use a finite‐difference time‐domain approximation to compute the acoustic Green's functions and perform a full waveform inversion for a multipole source, combining monopole and horizontal dipole terms. The inversion produces realistic estimates of VFR, on the order of 4 × 104 m3/s and well‐defined patterns of source directivity. This is the first application of acoustic waveform inversion at Mt. Etna. Our results demonstrate that acoustic waveform inversion is a mature and robust tool for assessment of source parameters and holds potential as a tool to provide rapid estimates of VFR in near real time.


Introduction
Over the past two decades infrasound has become increasingly popular as a tool for monitoring volcanic activity. A wealth of data collected at distances of few to thousands of kilometers from active volcanic vents has contributed to improving our ability to detect eruptions and track their evolution in real or near real time . Infrasound data gathered locally (within a few to tens of kilometers from of an active vent) are of interest in volcano monitoring for their potential to provide insights into eruption dynamics (Garcés et al., 1999;Johnson & Ripepe, 2011;Johnson et al., 2004;Kim et al., 2012;Matoza et al., 2014;Vergniolle & Ripepe, 2008). Recent studies have demonstrated the use of infrasound to evaluate parameters such as the time history of emissions during volcanic explosions (De Angelis et al., 2019;Johnson & Miller, 2014;Johnson et al., 2008;Kim et al., 2012Kim et al., , 2015, including calculations of volumetric and mass eruption rates, which are key parameters for evaluating the strength of eruptions. Infrasound is particularly attractive for volcano monitoring owing to its potential to provide estimates of eruption parameters even in adverse weather conditions, such as when cloud cover prevents measurements using other techniques and for its unmatched temporal resolution. These parameters are key inputs into models of atmospheric rise and transport of volcanic plumes; models linking infrasound-derived eruption source parameters to the rise dynamics and maximum height reached by volcanic plumes (Mastin et al., 2009;Woodhouse et al., 2013) have been tested with encouraging results (Caplan-Auerbach et al., 2010;Lamb et al., 2015;Ripepe et al., 2013;Vergniolle, 2008). It is worth noting that eruption parameters can also be recovered from Doppler radar measurements (Freret-Lorgeril et al., 2018), analysis of video recordings, satellite measurements, and gas-and-ash retrievals (Delle Donne et al., 2016;Dürig et al., 2015). The debate on the application of these techniques, including infrasound, and their associated uncertainty remains open .
Acoustic sources can be considered compact when their characteristic dimensions are significantly smaller than the wavelength of the radiated sound Kim et al., 2012); under this assumption, the acoustic wavefield produced by a point acoustic source can be represented as monopole (hemispherical radiation with equal strength in all directions), dipole (two monopoles of opposite phase), quadrupole (two dipoles), or a multipole (i.e., a combination of a monopole, dipole, and quadrupole; Lighthill, 1978;Rossing, 2015;Woulff & McGetchin, 1976). The importance of multipole sources is that they combine information on the magnitude of mass flux at the source during eruptions (monopole term) with estimates of the momentum imparted to the eruptive flow by external force fields such as those resulting from interaction with solid boundaries (dipole term). Dipole terms are commonly considered to provide estimates of source directivity, that is, how acoustic pressure varies with angular position with respect to the source. Early attempts to invert infrasound data for the assessment of eruption volume flow rate (VFR) adopted monopole sources, simple attenuation models, and data from individual sensors (e.g., Firstov & Kravchenko, 1996;Johnson & Miller, 2014;Johnson et al., 2008). Dipole sources were used to estimate gas velocity at Shishaldin (Vergniolle & Caplan-Auerbach, 2004), Augustine (Caplan-Auerbach et al., 2010), Eyjafjallajökull , and Redoubt (Lamb et al., 2015) Volcanoes. The first attempts to introduce a workflow for the simultaneous inversion of multistation infrasound data to resolve a multipole acoustic source were by Johnson et al. (2008) at Erebus (Antarctica) and Kim et al. (2012) at Tungurahua (Ecuador). Recent work by Kim and Lees (2014) and Lacanna and Ripepe (2013) demonstrated the influence of topography and vent geometry on the propagation of acoustic wavefields. Kim et al. (2015) developed a full waveform inversion technique, based on 3-D numerical Green's functions (GFs) computed by finite-difference time-domain (FDTD) modeling to estimate VFR from infrasound data. Fee et al. (2017) applied this technique for solving for a monopole source in order to estimate mass eruption rates at Sakurajima, Japan. More recently, Iezzi et al. (2019) extended the method, inverting infrasound waveforms for a multipole source, and tested its stability at Yasur Volcano, Vanuatu. They integrated data recorded by a microphone positioned above the active vents in order to investigate the 3-D sound radiation and its impact on the results of inversion. De Angelis et al. (2019) reviewed recent developments in volcano infrasound modeling and inversion, in particular, addressing linear acoustic wave theory and its application to the assessment of volcanic emissions using near-field acoustic data.
Here, we apply acoustic waveform inversion to infrasound data collected by a dense temporary deployment at Mt. Etna, Italy, during June-July 2017. Mt. Etna provides an ideal location to test and extend the applications of this method, due to its accessibility and nearly continuous eruptive activity. We compute 3-D numerical GFs based on the FDTD scheme of Kim and Lees (2014), to account for topographic effects, and perform full waveform inversions for a multipole source represented by the combination of a monopole and a horizontal dipole. We present estimates of VFR (i.e., monopole strength; Kim et al., 2015) and source directivity (i.e., dipole strength and orientation) and discuss the potential for application of acoustic waveform inversion in real time at Mt. Etna.

Data
We use data collected by a temporary infrasound network deployed at Mt. Etna, Italy, in the summer of 2017 ( Figure 1). The network consisted of 14 sensors (Chaparral M60-UHP), with a sensitivity of 9 mV/Pa and flat response between 0.03 and 245 Hz; data were sampled at 400 Hz with 24-bit resolution using DiGOS DATA-CUBE digitizers for a period of about 2 months. The network was centered around the summit area of Mt. Etna, with an azimuthal gap of about 50 • (Figure 1, and supporting information Figure S1). Stations were positioned at distances between 1 and 7 km from the active vent located within the New Southeast Crater ( Figure S1). During our campaign, volcanic activity remained relatively low. We recorded two distinct types of acoustic signals: (i) discrete explosions and (ii) tremor associated with background passive degassing. Here we focus on small-to-moderate size explosions recorded during 24-25 June and 8-9 July (supporting information Figure S2). During these periods the stations closer to the active vent recorded between 60 and 80 gas-and-ash explosions per day, with amplitudes rarely exceeding 5 Pa at 1 km from the source. Data from stations at greater distance from the vent exhibited low signal-to-noise ratio due to wind noise, making detection of small-amplitude signals challenging. The explosions analyzed were all gas-rich and produced small plumes reaching up to a few hundred meters above the vent, according to video footage and observer notes collected in the field. Additionally, daily SO 2 measurements for the duration of the experiment were provided by the Istituto Nazionale di Geofisica e Vulcanologia. This behavior is typical at Mt. Etna during periods of low-level activity when explosions are linked to bursting of isolated pockets of gas at shallow depth, instead of being driven by the influx of juvenile magma from depth (Pering et al., 2015). Explosion signals were characterized by a symmetric shape with an initial compression followed by a rarefaction phase and durations of between 2 and 5 s. The explosions radiated infrasound mainly in the 1.5-to 2.5-Hz frequency band, with most peaking at approximately 1.8 Hz (supporting information Figure S3). We identified a preliminary group of events suitable for further analysis by searching all stations on data filtered between 0.2 and 7 Hz, requiring that explosions were clearly recorded by at least five stations. A final suite of 14 high-quality explosions were selected for waveform inversion, including only signals with a signal-to-noise ratio (calculated as the ratio of the summed squares of 1 s of signal before and after the onset of an explosion) of at least 10.

FDTD GFs
The effects of topography on infrasound propagation over short distances have been extensively described in the literature (De Angelis et al., 2019;Kim & Lees, 2011;Lacanna et al., 2014;Matoza et al., 2009). In order to account for these effects, we calculated the atmospheric acoustic GFs using the 3-D FDTD numerical modeling scheme of Kim and Lees (2014) and Kim et al. (2015). The method splits the acoustic wave equation into two first-order differential equations for pressure and particle velocity and then approximates them using central differences in space and time with second-order accuracy. The perfectly matched layers technique sets absorbing boundaries to prevent outgoing waves from reflecting at the computational domain's boundaries. Topography was retrieved from a 2-m digital elevation model (DEM) from surveys conducted in 2017 over a 25-km 2 area around the summit craters, merged with a regional 30-m Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEM for the rest of the study area; the resulting ∼14.1 × 11.7-km area is parsed into a grid of 6 × 6-m cells. The DEM used was up-to-date with recent eruptive activity at the time of our deployment although, due to persistent cloud cover, some of the details of crater morphology remained unresolved. The FDTD method propagates the acoustic wavefield radiated by a multipole source (monopole and a horizontal dipole) Etna. Green's functions are plotted at increasing offset distance from the source (bottom to top; note that offset is not to scale). Green's functions represent the modeled acoustic pressure time series at different stations produced by propagating the acoustic wavefield generated by a source function (inset) over topography. (b) Modeled sound pressure level distribution measured in decibels relative to the reference atmospheric pressure of 70,108.54 Pa at 3,000 m above sea level. Topographical contour lines are plotted every 100 m. over the topography. Monopole and dipole GFs can be retrieved from a single simulation of a monopole source as the GF for an arbitrary dipole can be written in terms of the spatial derivative of the monopole's GF (Iezzi et al., 2019;Pierce, 1989).
The source impulse is approximated by a Blackman-Harris function with a cutoff frequency of 4.5 Hz (Figure 2). This frequency was selected to be sufficiently higher than the dominant frequency of all analyzed explosions (1.5-2.5 Hz) while taking into account the limitations imposed by the large computational domain grid. The source was positioned within the New Southeast Crater at the topographical surface; the position of the active vent was visible during the experiment and further confirmed by locations of explosion infrasound for all the events considered (supporting information Figure S4) using the semblance method, a grid search technique based on multichannel signal coherency (Almendros & Chouet, 2003;Cannata et al., 2011). The acoustic wavefield was propagated along the discretized topography using at least 20 grid points per wavelength, following recommendations of a 10-node grid per wavelength minimum (Wang, 1996). We defined a stratified atmosphere (Table S1 of the supporting information) based on density and temperature data recovered on 23 June 2017 by the nearest available radiosonde located in Trapani at 200 km from Mt. Etna (http://weather.uwyo.edu/ upperair/sounding.html). The GFs calculated for our network are shown in Figure 2a along with the spatial distribution of sound pressure level (the excess infrasound pressure relative to a reference atmospheric pressure of 70,108.54 Pa at 3,000 m above sea level) of the simulated acoustic wavefield (Figure 2b).

Waveform Inversion
Acoustic waveform inversion has been traditionally used to retrieve quantitative source models for volcanic explosions (De Angelis et al., 2019;Johnson et al., 2008;Kim et al., 2012). In this study we carried out a full waveform inversion following modifications to the method of Kim et al. (2015) and Fee et al. (2017), similar to Iezzi et al. (2019), which exploits 3-D GFs to invert for a multipole (monopole and dipole) source. The pressure generated by the combination of a monopole and a horizontal dipole can be written as (Pierce, 1989) follows: where p(x r , t) is the acoustic pressure recorded by a receiver located at position x r at time t, generated by a source located at x s . The * symbol indicates convolution; S(t) is the atmospheric mass flow rate; G m (x r , t; x s ) is the GF excited by a monopole source; are the GFs for a horizontal dipole; F x (t), F y (t) are the time histories of forces in the x and y directions. Equation (1) can then be rewritten in matrix form as follows: where d are the observed excess pressures (infrasound data) across the network, G is the matrix of the GFs, and m is the multipole source model. Equation (2) can be iteratively solved in a least square sense for m by minimizing the waveform misfit R:

Results
The explosions analyzed have similar waveforms characterized by compressional onsets followed by a rarefaction phase, and a short-duration coda ( Figure 3a); signals have overall durations of ∼4 s, and peak pressures between 0.5 and 7.5 Pa at the station closest to the active vent, ET04 (1.3 km from vent).
We performed waveform inversion for a multipole source on 14 selected signals to estimate the source time function (STF) of atmospheric mass flow rate (i.e., monopole strength) and source directivity (i.e., strength and direction of a horizontal dipole) over the STF. Prior to the inversion, waveforms were high-pass filtered with a cutoff frequency of 0.7 Hz, and data before the onset of the explosion were set to 0. In the results presented here, the STFs of atmospheric mass flow rates were transformed to atmospheric VFR using a standard atmospheric density of 0.909 kg/m 3 at 3,000 m above sea level. In order to reduce the effect of long-period drift that may influence source flow estimates, we follow the methods of Iezzi et al. (2019) and apply a postinversion linear detrend of the mass flow rate solution. The VFR for the 14 events analyzed are shown in Figure 3b. An example of application of the FWZPZF technique (Johnson & Miller, 2014) is illustrated in Figure 3c and compared with linear detrending (Iezzi et al., 2019). Details of the results of waveform inversion for one explosion are illustrated in Figures 3c and 3d, which show monopole and dipole strengths and the variation of dipole azimuth over the STF, respectively. The orientation of the horizontal dipole is stable (implications of this stability are further explained in the discussion), in the NW-SE direction, over the STF. For the 14 inversions, we obtained an average root mean square residual, R, of 0.17. The number of stations used varied between 10 and 14; for some events we discarded data from stations with signal-to-noise lower than 10 in order to avoid fitting noise during the inversion. An additional example of inversion results is provided in Figure 4, also showing the fit of the inversion model to the observed waveform data.
The range of peak VFR for the 14 explosions analyzed varied between and 18,000 and 96,000 m 3 /s. The cumulative volume produced by each explosion was obtained by integrating the VFR functions over the duration of the event. Cumulative volumes for the explosions in our data set ranged between 2.3 × 10 4 and 9.5 × 10 4 m 3 , which are of the same order of explosions with similar infrasound amplitudes recorded at Stromboli (Delle Donne et

Discussion and Conclusions
The use of infrasound waveform inversion methods for assessing eruption source parameters has become increasingly popular over the past decade (Fee et al., 2017;Johnson & Miller, 2014). Many previous studies have targeted large-amplitude, high-signal-to-noise ratio, Strombolian and Vulcanian style explosions, with excess pressures ranging between a few tens and few hundreds of pascals over relatively short distances from the source. Here we have successfully analyzed, for the first time, small-amplitude signals, thus extending the application range of this method. The use of 3-D numerical atmospheric acoustic GFs in our inversion accounts for the effects of topography on the propagation of the acoustic wavefield, addressing issues that may arise from the simplistic assumption that infrasound amplitudes decay according to geometrical spreading. Figure 2b clearly demonstrates the effects of topography on the acoustic wavefield, including significant signal attenuation in the NW sector due to the location of the active vent (within the New Southeast Crater) behind the topographic barrier represented by the summit North East, Voragine, and Bocca Nuova Craters (Figure 1).
Infrasound waveform inversion using numerical GFs, accounting for topographic effects, was performed in the past at Sakurajima Volcano (Japan), assuming a monopole source (Kim et al., 2015;Fee et al., 2017) and recently at Yasur Volcano, for a multipole source (Iezzi et al., 2019). Vulcanian activity at Sakurajima exhibits peak pressures of 10-400 Pa at a distance of 2.3 km from the vent and produce VFR of 6 × 10 5 -1.5 × 10 6 m 3 /s, 2 orders of magnitude larger than the small degassing explosions we observe at Etna. At Yasur, estimated VFR values also exceed our results, ranging from 1 × 10 5 to 6 × 10 5 m 3 /s. Values of volumetric flow rates comparable to those of our study were obtained for Strombolian activity at Erebus (1-2.5 × 10 4 m 3 ; Johnson et al., 2008) (Kim et al., 2012) were 1 order of magnitude larger than Etna. These studies, however, did not remove topographic effects from the inversion.
Our inversions produced very similar and stable dipoles oriented in the NNE-SSW direction (Figures 3 and  4) across the entire data set. The excellent azimuthal coverage and density of our network makes it well suited to resolve source mechanisms and acoustic radiation patterns. The similarity of source radiation patterns across the data set was expected owing to the repeating character of the selected explosion waveforms. However, the striking stability in dipole azimuth over the duration of the STF is uncommon in infrasound inversion studies (Iezzi et al., 2019), which tend to display a more erratic behavior reflecting variability in how sound is radiated at the source. We propose that the stability observed here may partly result from unresolved near-vent topography (Kim et al., 2012). We suggest that owing to insufficient resolution of the DEM in the crater area, some near-vent topography effects might not be fully accounted for in the numerical GFs.
At Etna the magnitude of the dipole source component appears, at least for the explosions considered in this manuscript, relatively small. To confirm this, we performed monopole-only inversions and compared the results to those of multipole inversion for the same events; differences in the overall shape and peak amplitudes of VFR functions were minimal. The difference in both peak VFR and cumulative erupted volume between the two workflows is of the order of 5% (supporting information Figure S5) with the monopole-only inversions providing slightly higher values. These observations suggest that, for the data considered here, a monopole source provides a realistic and sufficiently accurate estimate of VFRs (supporting information Figure S5). This is consistent with previous results from Iezzi et al. (2019) at Yasur and Fee et al. (2017) at Sakurajima. One important implication of this observation, if confirmed to have more general validity, is that if waveform inversion was implemented for use in real-time volcano monitoring, a monopole-only source would provide realistic estimates of atmospheric VFR.
Our results provide realistic time histories and peak values of VFR for the selected explosions in June-July 2017. The values obtained are consistent with information from visual observations and ground-based gas flux measurements, which are routinely performed at Mt. Etna (Pering et al., 2015;Tamburello et al., 2013 Although our results are encouraging, some limitations remain in the application of acoustic waveform inversion, in particular, in the assessment of the dipole source term. The vertical component of the radiated sound field and its effect on the assessment of dipole strength appears only when sensors are located over the source (i.e., tethered to balloons over the crater; Iezzi et al., 2019), and effects of unresolved topography may overprint true source directivity. Furthermore, we note the challenge of translating atmospheric VFR into mass eruption rates (Fee et al., 2017), a more useful parameter for input into models of atmospheric ash transport and dispersal. Conversion from volume to mass flow rate relies on measurements of the density of eruptive flows that, in turn, depend on detailed knowledge of the composition of erupted mixtures, which is the relative proportions of ash and gas in eruptive plumes . Such measurements are still comparatively rare at volcanoes, and, when available, the temporal resolution of such data is still considerably lower than infrasound. At Mt. Etna, Doppler radar measurements have been recently used to estimate mass eruption rates (Freret-Lorgeril et al., 2018) and hold promise for future use alongside infrasound and other multiparameter data streams. Other potential approaches to estimate flow density and gas composition include video analysis of the transition between the gas-thrust and buoyancy phases (Dürig et al., 2015) and gas measurements (Dalton et al., 2010;Delle Donne et al., 2016;Tamburello et al., 2013).
In this study we have demonstrated the feasibility of waveform inversion to retrieve time histories of atmospheric VFR for weak explosive activity at Mt. Etna. Our results, alongside an extensive body of literature on infrasound inversion, suggest that the method is mature for real-time implementation and integration within the volcano monitoring program at Mt. Etna and other volcanoes where infrasound data are available.
Owing to the relatively high computational cost, we envision that the implementation of infrasound inversion in real or near real time should be based on precomputed libraries of acoustic GFs, similar to routine practices in seismic moment tensor inversions. We surmise that acoustic waveform inversion holds potential to help volcano monitoring owing to its ability to provide real-time estimates of eruption source parameters, which can be used to assess the magnitude of explosions and as inputs into models of atmospheric rise of volcanic ash plumes for larger explosive activity.