Uncertainty in Detection of Volcanic Activity Using Infrasound Arrays: Examples From Mt. Etna, Italy

ORIGINAL RESEARCH published: 28 May 2020 doi: 10.3389/feart.2020.00169 Frontiers in Earth Science | www.frontiersin.org 1May 2020 | Volume 8 | Article 169 Edited by: Andrea Cannata, University of Catania, Italy Reviewed by: Stephen Arrowsmith, Southern Methodist University, United States Stephen R. McNutt, University of South Florida, United States *Correspondence: Silvio De Angelis [email protected] Specialty section: This article was submitted to Volcanology, a section of the journal Frontiers in Earth Science Received: 13 February 2020 Accepted: 01 May 2020 Published: 28 May 2020 Citation: De Angelis S, Haney MM, Lyons JJ, Wech A, Fee D, Diaz-Moreno A and Zuccarello L (2020) Uncertainty in Detection of Volcanic Activity Using Infrasound Arrays: Examples From Mt. Etna, Italy. Front. Earth Sci. 8:169. doi: 10.3389/feart.2020.00169 Uncertainty in Detection of Volcanic Activity Using Infrasound Arrays: Examples From Mt. Etna, Italy Silvio De Angelis1*, Matthew M. Haney2, John J. Lyons2, Aaron Wech2, David Fee 3, Alejandro Diaz-Moreno4and Luciano Zuccarello5,6 1School of Environmental Sciences, University of Liverpool, Liverpool, United Kingdom, 2Alaska Volcano Observatory, U.S. Geological Survey Volcano Science Center, Anchorage, AK, United States, 3Alaska Volcano Observatory, Geophysical Institute, University of Alaska Fairbanks, Fairbanks, AK, United States, 4School of Geography, Politics and Sociology, Newcastle University, Newcastle upon Tyne, United Kingdom, 5Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Pisa, Pisa, Italy, 6Departamento de Fisica Teorica y del Cosmos, University of Granada, Granada, Spain The injection of gas and pyroclastic material from volcanic vents into the atmosphere is a prolific source of acoustic waves. Infrasound arrays offer efficient, cost-effective, and near real-time solutions to track the rate and intensity of surface activity at volcanoes. Here, we present a simple framework for the analysis of acoustic array data, based on least-squares beamforming, that allows to evaluate the direction and speed of propagation of acoustic waves between source and array. The algorithms include a new and computationally efficient approach for quantitative assessment of the uncertainty on array measurements based on error propagation theory. We apply the algorithms to new data collected by two 6-element infrasound arrays deployed at Mt. Etna during the period July–August 2019. Our results demonstrate that the use of two infrasound arrays allowed detecting and tracking acoustic sources from multiple craters and active vents associated with degassing and ash-rich explosions, vigorous and frequent Strombolian activity, opening of new eruptive fractures and emplacement of lava flows. Finally, we discuss the potential use of metrics based on infrasound array analyses to inform eruption monitoring operations and early warning at volcanoes characterized by episodic intensification of activity. Keywords: volcano infrasound, infrasound arrays, volcanic degassing, ash explosions, lava flow, Mt. Etna 1. INTRODUCTION Over the past two decades infrasound has become an increasingly popular tool to monitor volcanoes (e.g., Fee and Matoza, 2013; McNutt et al., 2015). Because of its ability to detect and discriminate shallow and sub-aerial volcanic activity (Matoza et al., 2019) infrasound is a desirable complement to seismology for monitoring unrest, and to detect and track the evolution of eruptions in realor near real-time (e.g., Ripepe et al., 2009, 2018; Cannata et al., 2013; Coombs et al., 2018) that is over time scales of the order of few seconds to few minutes. Acoustic waves are generated when the atmosphere is perturbed from equilibrium (e.g., Garcés et al., 2009). At volcanoes, small and large explosions, gas-and-ash jets and plumes, sector collapses, rockfalls, pyroclastic flows, and lahars, are likely to generate infrasound over a wide range of magnitudes, frequencies, and durations (Johnson and Ripepe, 2011; Fee and Matoza, 2013; Johnson and Palma, 2015; McNutt et al., 2015; Matoza and Fee, 2018). Infrasound is recorded by band-sensitive microphones at different scales De Angelis et al. Volcano Infrasound Array Analysis (Fee and Matoza, 2013) from local (<20 km) to regional (several 10s to few 100s of km), and even global (several 100s to 1,000s of km). At the local scale, microphones are deployed either as individual sensors within distributed networks or as smallaperture arrays (i.e., tight clusters of three or more instruments at 50–150 m from one another) at distances of between few hundreds of meters to within 20 km of an active volcanic vent. Data from distributed networks have traditionally been used for absolute source location and event discrimination (e.g., Cannata et al., 2009), and to evaluate eruption source parameters (e.g., Caplan-Auerbach et al., 2010; Fee et al., 2017; De Angelis et al., 2019; Diaz-Moreno et al., 2019). On the other hand, the close spacing between sensors deployed as a smallaperture array allows detection of coherent infrasound even in low signal-to-noise ratio (SNR) conditions, characterization of the direction-of-arrival (DOA), and the apparent speed of propagation of acoustic waves. Notably, the use of infrasound arrays even at distances of a few tens of kilometers from active vents provides robust and efficient remote detection of eruptive activity and source discrimination (e.g., Fee et al., 2010), thus reducing risks for observatory personnel during field campaigns. Johnson (2004) and Ripepe et al. (2007) discussed applications of infrasound array methods to track degassing and eruptive activity from multiple active vents at Stromboli volcano (Italy). Ripepe et al. (2018) reported on over a decade of acoustic monitoring at Mt. Etna (Italy) with two small aperture infrasound arrays, demonstrating the ability to detect eruption onset in real-time, discriminate source position, and dispatch rapid notifications of the ongoing activity to local civil protection authorities. Matoza et al. (2007) analyzed array data during activity at Mt. St. Helens (USA) from October 2004 to March 2005; they discussed the key role of infrasound in separating surface from deeper processes and in identifying the timing and assessing the magnitude of eruptive events. Previous studies focussing on (local scale) infrasound array applications to volcano monitoring do not provide quantitative estimates of the true uncertainty associated with their measurements. In this manuscript we tackle this issue presenting a simple framework for infrasound array processing based on least-squares beamforming (Olson and Szuberla, 2005; Haney et al., 2018) and introducing a new scheme for quantitative assessment of the uncertainty on estimates of DOA and apparent velocity. We test the method on data collected by two arrays deployed at Mt. Etna (Italy) during the summer of 2019 (Figure 1). We describe the array processing workflow, and show how infrasound data can be used to discriminate and track volcanic activity, from background degassing to individual ashrich explosions, persistent Strombolian activity and lava effusion. 2. ACTIVITY AT MT. ETNA: JULY–AUGUST 2019 Mt. Etna, Italy, is one of the most active volcanoes in the world. Several eruptive episodes have taken place in the past two decades, at its summit, from the North East (NEC), Voragine (VOR), Bocca Nuova (BN), and New South East (NSEC) craters. Since 2011, more than 50 effusive events have taken place in the area, which is visited by thousands of tourists every year (e.g., De Beni et al., 2019; Sciotto et al., 2019). During the last decade lava fountaining has frequently been observed from the NSEC, as well as shorter episodes of Strombolian activity from BN, VOR, and the NEC (Sciotto et al., 2019). Eruption at Etna occasionally evolves into episodes of more intense activity referred to as paroxysms that is, a significant increment in the rate and intensity of explosions from one or multiple active vents accompanied by emplacement of lava flows and/or generation of significant ash plumes. Activity at Mt. Etna during July– August 2019 was marked by the occurrence of two paroxysms on 18 and 27–28 July, accompanied by the emplacement of lava flows from lateral vents in the NSEC area (Figure 1). The weeks preceding each paroxysm were characterized by both sustained and explosive degassing occasionally punctuated by ash-rich explosions (Figures 2A,B). The Istituto Nazionale di Geofisica e Vulcanologia (INGV) reported intense degassing and observations of four large explosions from the NEC on 2, 3, 8, and 13 July (Figure 2B). Degassing activity from VOR was low-to-moderate, whilst BN produced deep intra-crater gas explosions with minor amounts of ash (Figure 2A). Activity at the NSEC increased gradually from vigorous degassing to ash-rich Strombolian activity throughout the first week of July; degassing-only activity resumed on 7 July and escalated again into strong and nearly continuous Strombolian explosions on 18 July. Explosion rates varied from one every 1–2 min to several per minute, eventually leading to the opening of a fracture and subsequent emplacement of a lava flow on the NE flank of the NSEC (23:09 UTC on 18 July). During the night on 19 July, Strombolian activity shifted from NSEC to NEC, and gradually decreased until it completely halted—including lava flow effusion—in the morning of 20 July. On 25 July, a new phase of eruption began at the NSEC with the onset of Strombolian activity, transitioning into nearly continuous lava fountaining early in the morning of 27 July. This activity was accompanied by sustained ash emissions (Figures 2C,D) just a few hours before the opening of two new vents on the southern flank of the NSEC from which two new lava flows developed (Figure 1). A strong explosion occurred on July 27 at 12:21 UTC, accompanied by an ash plume that reached nearly 4 km above the vent (Figure 2C). Strombolian activity ceased on 28 July at 03:40 UTC, while lava flows continued until late that day. The remainder of our deployment period was characterized by background levels of degassing from NSEC, VOR, and BN, while NEC continued producing episodic ash explosions during August. INGV reported no major changes in deformation during the 2 months of our temporary deployment, while daily gas emissions increased notably during 15–29 July and then returned to background levels. Seismic tremor fluctuated with marked peaks on 6 July, coincident with the increase in Strombolian activity at NSEC. Seismic tremor remained stable, at high levels, during both paroxysms, and eventually returned to background during August. Finally, infrasound locations provided by INGV successfully identified the first paroxysm on July 18, but high Frontiers in Earth Science | www.frontiersin.org 2May 2020 | Volume 8 | Article 169 De Angelis et al. Volcano Infrasound Array Analysis FIGURE 1 | (A) Map showing: (1) the locations of two infrasound arrays (green triangles), ENEA and ENCR, deployed at Mt. Etna, Italy, during July–August 2019; (2) approximate locations of vents active in the summit area (white circles of variable size according to vent dimensions as qualitatively inferred from visual observations) during the deployment period; (3) active lava flows during the deployment period. (B,C) Details of the geometry of the two infrasound arrays, ENEA and ENCR. FIGURE 2 | Images of eruptive activity at Mt. Etna, Italy, during July–August 2019: (A) Intense sustained degassing activity from the North East Crater (NEC) and deep intra-crater gas-rich explosions from Bocca Nuova (BN) recorded during a UAV flight on 5 July 2019 at 08:00 UTC (Photo credit: A. Diaz-Moreno, University of Liverpool); (B) Ash-rich explosion from NEC on 3 July 2019 at 10:06 UTC (Photo credit: A. Diaz-Moreno, University of Liverpool); (C) Explosive activity accompanied by sustained ash emissions observed on 27 July 2019 at 12:21 UTC from the New South East crater (NSEC) area (Photo credit: Giuseppe Di Stefano/Marco Restivo, Etna Walk); (D) Intense Strombolian activity and lava flow originating from a vent in the NSEC area on 27 July 2019 (Photo credit: INGV-Osservatorio Etneo). Frontiers in Earth Science | www.frontiersin.org 3May 2020 | Volume 8 | Article 169 De Angelis et al. Volcano Infrasound Array Analysis wind conditions on 27–28 July prevented the network from detecting the second paroxysm. The activity recorded during the summer of 2019 is typical at Mt. Etna; other eruptions with a similar fingerprint have been reported in the past (Corsaro et al., 2017; De Beni et al., 2019; Polacci et al., 2019). 3. DATA In this study, we used data recorded by two 6-element infrasound arrays (ENEA and ENCR, Figure 1), deployed at Mt. Etna between 2 July and 25 August 2019. Changes in eruptive activity during this period, summarized in section 2 of this manuscript, are reported in the INGV activity bulletins and were confirmed by both Unmanned Aerial Vehicle (UAV) imagery (Figure 2A) and ground-based visual observations (Figures 2B–D) gathered during the deployment by the authors. The ENEA and ENCR arrays were installed at distances of between 1 and 1.8 km from the active vents (Figure 1). Each had similar configurations, deployed with an ∼100 m aperture (i.e., the largest distance between any two elements within an array) pentagon shape with a central element. Care was taken, considering the constraints imposed by topography, that the difference in elevation between any two microphones within each array was small; for both arrays, this difference did not exceed ∼30 m (see Methods section for additional details). The arrays were designed to have apertures large enough to discriminate acoustic arrivals and small enough to record coherent infrasound across all microphones (Ripepe et al., 2007). Data from both arrays were recorded onsite with a sampling frequency of 100 Hz at 24-bit resolution using DiGOS Data-cube digitizers (https://digos.eu/seismologyand-cubes/); ENEA was equipped with Chaparral Model 60 Ultra High Pressure microphones (full-scale range of 2000 Pa peak-topeak, flat response between 0.03 and 245 Hz), and ENCR with IST2018 sensors (full-scale range of 480 Pa peak-to-peak, flat response between 0.1 and 40 Hz, as described in Grangeon and Lesage, 2019). Pressure amplitudes recorded by the two arrays ranged from few Pa for signals associated with both sustained and explosive degassing, to several tens of Pa during intense and persistent Strombolian activity. Examples of multi-channel data recorded by both ENEA and ENCR are shown in Figures 3A–E; these include an ash-rich explosion from NEC (Figure 3C) and repeated deep intra-crater explosions from BN (Figure 3D) on 2 July 2019, as well as Strombolian activity from NSEC on 18 July 2019 (Figure 3E). It is worth noting how infrasound waveforms in Figures 3B,E exhibit a marked asymmetry, although the investigation of this intriguing feature falls beyond the scopes of this manuscript. 4. METHODS Here we discuss an array processing workflow to derive estimates of DOA and horizontal velocity from volcano infrasound array data. The method solves, via least-squares inversion, the problem of fitting a plane wave arrival traveling from an azimuth θwith horizontal velocity vto a vector of measured time delays between traces across the array (e.g., Claerbout, 1986; Olson and Szuberla, 2005; Haney et al., 2018). We also introduce a new method to quantify errors on the measurements of θand vusing standard propagation of error. For an array of nsensors, a data vector d=(δt1,δt2,...,δtN) of N=n(n−1)/2 delay times between the elements of the array can be estimated, for instance using waveform cross-correlation; a linear relationship exists between dand a model vector, m=(sx,sy)T=(sin θ/v, cos θ/v)Tof slowness (defined as the inverse of velocity) in the East-West (sx) and North-South (sy) directions: d=Gm (1) where Gis a Nx2 matrix of horizontal distances between all pairs of array elements. When the measurements of delay times are affected by error, Equation (1) can be re-written to explicitly include it as: d=Gm +ǫ(2) where ǫis the vector of errors in the estimates of time delays. These errors can be assumed to be normally distributed, with zero mean and variance σ2 τ(Olson and Szuberla, 2005). The solution to (2) for mis then found by minimizing the sum of squared errors: E=ǫǫT=(d−Gm)T(d−Gm) (3) that is: m=(GTG)−1GTd(4) Finally, DOA and apparent trace velocity across the array can be estimated from the solution vector mas: (v=1/qs2 x+s2 y θ=tan−1(sx/sy)(5) It is important to note that the beamforming analysis in (1– 5) is carried out in the horizontal plane (Edwards and Green, 2012), that is under the assumption that the contribution of differences in sensor elevations to the time delay measurements can be neglected; for this reason it is crucial that such differences are small. For the most part, studies that apply array processing to volcano infrasound data either provide qualitative statements on the uncertainty associated with estimates of DOA and trace velocity or discuss the theoretical azimuthal resolution of the array in terms of its aperture in relation to the wavelength of the signals analyzed (e.g., Ripepe et al., 2018). The statistical confidence in the estimates of θand vare rarely discussed (e.g., Szuberla and Olson, 2004); precise estimates of DOA and trace velocity are made difficult by the ubiquitous presence of noise in the data and the complex propagation of the acoustic wavefield, which affect measurements of time delays across arrays. Under the assumption that errors on the measurements of time delays across an array are normally distributed, with zero mean and variance σ2 τ, the deviation of mfrom its expected value, due to Frontiers in Earth Science | www.frontiersin.org 4May 2020 | Volume 8 | Article 169 De Angelis et al. Volcano Infrasound Array Analysis FIGURE 3 | Infrasound data recorded by the ENEA (blue) and ENCR (black) arrays during July 2019: (A) 5 min of continuous infrasound activity recorded by ENEA on 2 July 2019; (B) 6 h of intense Strombolian activity recorded by ENCR on 18 July 2019; (C,D) Detail (gray shaded areas in A) of ENEA infrasound associated with an ash-rich explosion from the NEC and deep intra-crater gas explosions from BN, respectively; (E) Detail (gray shaded area in B) of ENCR infrasound recorded during vigorous Strombolian activity on 18 July 2019. Note the difference in infrasound amplitudes between signals in (A,B). errors in the estimates of time delays, is described by the model covariance matrix, C(m), given by: C(m)=σ2 τ(GTG)−1=σ2 τD(6) The matrix C(m) is symmetrical and can be written in terms of the variances (σ2 sx,σ2 sy) and covariance (σ2 sx,sy) of sxand sy: C(m)="σ2 sxσ2 sx,sy σ2 sx,syσ2 sy#(7) Frontiers in Earth Science | www.frontiersin.org 5May 2020 | Volume 8 | Article 169 De Angelis et al. Volcano Infrasound Array Analysis FIGURE 4 | Array processing of 35 min of data recorded at ENEA on 2 July 2019 between 9:50 and 10:25 UTC. Note (bottom panel) consistent detection of deep intra-crater activity from BN briefly interrupted by a larger ash-rich explosion from NEC (see section 2). Equations (6) and (7) provide a link between the variances and covariance of sxand sy, and the variance of time delay measurements, that is:      σ2 sx=σ2 τD(1, 1) σ2 sy=σ2 τD(2, 2) σ2 sx,sy=σ2 τD(1, 2) (8) We are interested in error estimates on DOA and trace velocity, which are non-linear functions of sxand syas shown in (5). In order to estimate the variances for θand v(σ2 θand σ2 v, respectively) we apply standard propagation of errors theory (e.g., Vardeman and Jobe, 2001) to (5). Neglecting high-order terms we obtain:    σ2 v≈σ2 sx∂v ∂sx2 +σ2 sy∂v ∂sy2 +2σ2 sx,sy∂v ∂sx∂v ∂sy σ2 θ≈σ2 sx∂θ ∂sx2 +σ2 sy∂θ ∂sy2 +2σ2 sx,sy∂θ ∂sx∂θ ∂sy(9) The terms in (9) are obtained by simply differentiating (5) and evaluating them at the values for sxand sygiven by the solution of (4), that is: (σ2 v≈σ2 sxs2 xv6+σ2 sys2 yv6+2σ2 sx,sysxsyv6 σ2 θ≈σ2 sxs2 yv4+σ2 sys2 xv4+2σ2 sx,sysxsyv4(10) Finally, the standard deviations of trace velocity and DOA, σvand σθ, are obtained by taking the square root of (10). 5. RESULTS We applied the array processing workflow described in section 4 of this manuscript to continuous infrasound data recorded at Mt. Etna between 2 July 2019 and 25 August 2019. Data collected by two infrasound arrays, ENEA and ENCR, were pre-processed by applying a band-pass (2-pole, zero-phase) filter within the frequency band of interest, which for the activity recorded during our experiment was between 0.7 and 15 Hz. The delay times between pairs of sensors across each array were determined using cross-correlation with sub-sample accuracy (e.g., Haney et al., 2018) within 10-s sliding windows overlapping by 5 s. DOA and trace velocity were calculated for all data windows by leastsquares inversions according to (1–5), and their variances were calculated according to (7–9). The mean of the normalized crosscorrelation maxima (MCCM) within each data window reflects the level of coherence between signals across the array; only values of DOA and trace velocity corresponding to MCCM > 0.5 were considered to ensure that inversion was performed for coherent signals across the arrays. Examples of estimates of DOA and trace velocity, along with their uncertainties, are shown in Figures 4,5. In Figure 4 we Frontiers in Earth Science | www.frontiersin.org 6May 2020 | Volume 8 | Article 169 De Angelis et al. Volcano Infrasound Array Analysis FIGURE 5 | Array processing of 51 h of data recorded at ENEA and ENCR, starting on 18 July 2019, 00:00 UTC: (A) ENCR shows stable detections during Strombolian activity (18:00–24:00 on 18 July and 12:00–22:00 on 19 July) as well as during a period of lava effusion from the NE flank of NSEC (23:09 on 18 July, see section 2 in this manuscript); (B) ENEA captures the most intense Strombolian activity from the NSEC but shows larger fluctuations in DOA (direction of arrival) than ENCR. This is due to less favorable signal-to-noise ratio for activity at the NSEC as well as its proximity to secondary sources within the BN crater. Both arrays detect a clear shift in activity from NSEC to NEC late during the night between 19 and 20 July. Note how uncertainties and scatter in locations are consistently smaller at the array closer to the active vent. show data recorded at ENEA on 2 July 2019 (the only array installed at this time). Figure 5 illustrates results from both ENEA and ENCR for just over 2 days of activity starting on 18 July 2019. These examples are representative of the variable activity, during July–August 2019, from multiple vents (Figure 1) and across all summit craters at Mt. Etna. Figure 4 shows stable Frontiers in Earth Science | www.frontiersin.org 7May 2020 | Volume 8 | Article 169 De Angelis et al. Volcano Infrasound Array Analysis locations in the direction of BN associated with deep intracrater explosions (Figures 2A,3D) and a sharp change in DOA at about 10:06 UTC, which corresponds to a larger ash-rich explosion from the NEC (Figures 2B,3C); trace velocities are of the order of 330–335 m/s with standard deviations of 5–15 m/s, and uncertainties on DOA estimates of 0.5–2◦(≈10–40 m at the distance between ENEA and the vent). Locations in Figure 5 are dominated by activity at NSEC, in particular between 18:00 and 24:00 on 18 July and during 12:00–19:00 on 19 July, when frequent and intense Strombolian activity was observed at this crater (Figures 3B,E). Values of trace velocities at ENEA during this period are of the order of 310–335 m/s, with uncertainties of 5–15 m/s and standard deviations on DOA between 1 and 2◦ (≈30–60 m); at ENCR, trace velocities are in 330–335 m/s range with uncertainties of 2–8 m/s, and uncertainties on DOA of 0.5– 2◦(≈10–40 m). Notably, both arrays detected a shift in activity in the direction of the NEC, also reported by INGV at about 23:00 on 19 July (section 2). The two arrays were installed to provide the best coverage of all active craters at Mt. Etna (within limitations imposed by site access and safety). Owing to their positions relative to the active vents ENEA was ideally located to independently discriminate activity from NEC, whereas ENCR was positioned to optimize detection of activity from the NSEC area (Figure 1). In Figure S1, we have provided an additional example of array locations showing the stability and low uncertainties associated with detections at ENCR during a period of Strombolian activity, lava fountaining and sustained ash emissions at the NSEC on 27–28 July 2019 (see section 2 in this manuscript, Figures 2C,D). 6. DISCUSSION AND CONCLUSIONS We have presented an algorithm for the inversion of infrasound array data that allows estimating DOA with apparent sound velocity, and their associated uncertainties. We have applied the workflow to data collected by two 6-element infrasound arrays deployed at Mt. Etna during the summer of 2019; continuous detections from the two arrays tracked the activity observed at the volcano, including shifts in degassing and eruptions across multiple summit craters and vents. The results presented in Figures 4,5, suggest that: (i) array location relative to a complex system of active vents is key to allow discrimination of variable activity across multiple craters, and (ii) the quality of detection and thus the final estimates of uncertainties are crucially influenced by SNR levels. Figure 5, for instance, shows higher quality (i.e., lower uncertainties) detections from ENEA for activity located at the NEC, while ENCR seems better suited to track eruptions in the NSEC area; ENEA is also able, albeit with comparatively large uncertainties, to track small intracrater explosions from BN (Figure 5B) even during periods of elevated Strombolian activity at the NSEC. At the first order, this observation arguably reflects variable SNR at the two sites, resulting from the interplay between array proximity to the source (NEC and BN closer to ENEA; NSEC closer to ENCR) and wind noise levels (qualitatively observed to be consistently higher at ENEA). Wind strength and direction are factors that can potentially introduce a bias in the estimates of DOA and trace velocity from infrasound arrays (e.g., Schwaiger et al., 2020). At large source-receiver distances time reversal location of infrasound sources, propagating the acoustic wavefield in a windy atmosphere, has been used to assess the misfit between known source backazimuths and DOA estimates (e.g., Lonzaga, 2016). At the local scale, Johnson et al. (2012) demonstrated how temporal variations in acoustic parameters, such as infrasound travel times over short distances—<20 km—can be exploited to infer atmospheric conditions, including the strength of wind in a vertically stratified atmosphere. In theory, appropriate corrections for the effect of wind on infrasound measurements could be introduced in array inversion workflows; the main challenge for this lies in the fact that wind measurements or models with the required temporal and spatial resolution for the local scale are generally not available. In addition to wind noise, uncertainty on estimates of DOA and sound velocity are further linked to array configuration and its position relative to the sources; these factors control the degree to which the measured time delays across the array correspond to a physically realizable set of delays associated with the propagation of a plane wave across the array (Szuberla and Olson, 2004). Qualitatively, the plane wave approximation is considered valid at sourceto-array distances much greater than the aperture of the array (e.g., Almendros et al., 1997), a condition met in our study for all source-array combinations. Finally, effects from topography, such as diffraction can affect acoustic propagation and introduce additional bias in array estimates of DOA and sound speed. At the local scale Fee et al. (2019) recently investigated the effects of topography on back-projection locations of volcano infrasound sources integrating improved travel time calculations from finite difference modeling of the acoustic wavefield. At Mt. Etna, specifically, Diaz-Moreno et al. (2019) performed 3D finite difference simulation of acoustic wavefield propagation; their results suggest that the effects of topography on acoustic propagation from sources in the area of the summit craters to the locations of the arrays used in this study are minor. A quantitative estimate of how well the measured time delays correspond to a plane wave crossing the array is directly provided by the a priori variance of the time delay measurements, which is given by (Szuberla and Olson, 2004): σ2 τ=R2 0 N−r(11) where (following the same notation of Equations 1–10) R2 0= dT(I−R)d, with Ibeing the identity matrix and R= G[(GTG)−1−1]GT; N is the number of station pairs and ris the rank of R. We calculated σ2 τfor periods with high MCMM >0.5; values were of the order of 10−2s, that is ≈1/fs(with fsbeing the data sampling frequency) as also reported in other studies (e.g., Szuberla and Olson, 2004). This gives confidence that the plane wave assumption does not introduce significant bias on the results of our array analysis. On the other hand for data with low MCMM, that is no coherent signal traveling across the array, the values of σ2 τwere very high, up to a few seconds. Frontiers in Earth Science | www.frontiersin.org 8May 2020 | Volume 8 | Article 169 De Angelis et al. Volcano Infrasound Array Analysis FIGURE 6 | One-minute Infrasound Parameter during July–August 2019 for arrays ENEA (A) and ENCR (B); periods when no data were available are also indicated. Detail (zoom in of gray areas in A,B) of Infrasound Parameter values at ENEA (C) and ENCR (D) for the period 17–21 July 2019. The figure shows an obvious increase in both occurrence (i.e., increasing array detections per minute) and intensity (increasing signal amplitude) of Strombolian activity; the Infrasound Parameter peaks in coincidence with the opening of a new fracture and emplacement of a lava flow on the NE flank of the NSEC, late during the night on 18 July (see sections 2 and 6). In our workflow, uncertainties on apparent velocity and DOA are evaluated using the theory of error propagation. While the variances estimated using this method are a first order approximation to the true uncertainty, this solution is easy to implement with very low computational overhead, an obvious advantage for deployment within real-time volcano monitoring systems. An exact solution for the statistical confidence in the estimates of DOA and signal velocity for planar arrays (i.e., without significant variations in elevation across the array) is discussed in Szuberla and Olson (2004). We benchmarked the results of our workflow against this solution and found good agreement between the two approaches (Figure S2). Our study also reveals that combining estimates of DOA from multiple arrays improves the ability to discriminate multiple active vents. The ENEA and ENCR arrays were ideally positioned to independently discriminate sources within the NEC and NSEC, respectively (Figure 1). In the absence of complementary observations, however, resolving activity from variable vents within the other two summit craters (BN and VOR) frequently requires joint measurements from two arrays; combining the DOAs obtained from both arrays, a source location could be uniquely identified (Figure S3). Finally, we note that our workflow implicitly assumes a single arrival within each signal window analyzed. At volcanoes with multiple vents separated by a small distance different sources could be simultaneously active, and thus, produce multiple arrivals within a signal analysis window. This is discussed, for example, by Yamakawa et al. (2018) at Stromboli; the authors demonstrate how the MUSIC algorithm can help to resolve multiple active sources. We find that when active vent are separated by comparatively large distances and have different characteristic dimensions, such as Mt. Etna, appropriate selection Frontiers in Earth Science | www.frontiersin.org 9May 2020 | Volume 8 | Article 169