Acoustic observations of lava fountain activity during the 2021 Fagradalsfjall eruption, Iceland

The 2021 eruption within the Fagradalsfjall volcanic system in Iceland provided a rare opportunity to record acoustic data generated by a basaltic fissure. Eruptive activity in May 2021 was defined by a sequence of repetitive lava fountaining activity. Here we describe key observations and analysis conducted on acoustic data recorded by a four-element infrasound microphone array near the eruption site. Detailed inspection of acoustic waveforms and comparisons with seismic data and lava fountain height measurements revealed a complex eruptive sequence during each lava fountain event: acoustic tremor during peak lava fountaining was followed by a transition to Strombolian-style activity with distinct high-amplitude impulsive waveforms. Quantitative comparisons to jet noise spectra find complex turbulence acoustics during each event, with evidence of variations in the wavefield centred on peak lava fountain heights. Strombolian explosions could mostly be modelled by oscillations of bursting gas slugs at the top of the magma column, with a minor number of events exhibiting Helmholtz resonance behaviour instead. We find an increase in bubble radii between early and late May, suggesting a widening of the upper conduit during the lava fountain sequence. Finally, we propose that higher acoustic amplitudes, in addition to a wider conduit in late May, indicate higher gas flux through the conduit culminating in shorter lava fountain events. This study highlights the value of deploying acoustic sensors for providing additional constraints on eruption dynamics and source parameters during effusive fissure eruptions in Iceland and elsewhere.


Introduction
Atmospheric acoustic signals with frequencies ranging from 0.02 to 20 Hz are classified as infrasound and are regularly documented during volcanic activity (Johnson and Ripepe 2011;Fee and Matoza 2013). These atmospheric perturbations may be generated by processes including shortduration explosions, sustained volcanic jets, and mass flows . Localisation and quantification of acoustics during eruptions can provide information on shallow processes within the conduit and above the vent. Therefore, the deployment of microphones around an active volcano can provide additional insights into activity not readily available by other means.
Relative to other types of eruptive activity such as Strombolian and Vulcanian, acoustic records of basaltic fissure eruptions with Hawaiian-style lava fountains are rare with previous studies limited to events at Etna (Gresta et al. 2004;Cannata et al. 2009;Cannata et al. 2011;Ulivieri et al. 2013;Ripepe et al. 2018) and Kīlauea volcanoes (Fee et al. 2011;Lyons et al. 2021;Gestrich et al. 2021). Acoustic activity during these events is usually manifested as semi-continuous tremor with broadband frequency spectra characteristics. Using sensor arrays and networks, it is possible to track the activity progression along fissures, including the opening of new vents (Cannata et al. 2011;Fee et al. 2011), and quantifying the effusion rate of high-speed, highly channelised lava flows (Lyons et al. 2021). Acoustic source processes at the vent(s) during fissure eruptions have been modelled as either Strombolian-style explosions (Cannata et al. 2011) or a sustained jet flow (Gestrich et al. 2021).
Low viscosity magma allows large gas bubbles to ascend through the magma column and manifest at the surface as distinct, impulsive, high-amplitude Strombolian bursts. During this type of activity, at least two different acoustic source processes have been described. Firstly, strong gas bubble oscillation immediately prior to bursting at the lava surface generates an m-shaped waveform where the positive peaks are less intense than the negative peak (Vergniolle and Brandeis 1994;Vergniolle and Ripepe 2008). Secondly, longer harmonic waveforms with an exponentially decaying coda are linked to Helmholtz resonance of gas within or above the bubble Cannata et al. 2009;Fee et al. 2010b;Goto and Johnson 2011). Modelling of Strombolian acoustics allows for quantitative estimates of parameters such as bubble radius, bubble length, and overpressure . During the 2008 flank fissure eruption at Etna volcano, both the bubble oscillation and Helmholtz resonance source model were detected (Cannata et al. 2009). For different vents across the fissure, waveform inversions for bubble oscillations produced estimated bubble radii of 2-3.5 m with bubble lengths of 7.2-7.6 m (Cannata et al. 2009;Cannata et al. 2011). In contrast, events generating Helmholtz resonance were modelled by bubbles of 6 m radius and up to 40 m length (Cannata et al. 2009).
Recent studies have demonstrated how eruptions can generate an infrasonic form of 'jet noise' such as those generated by small-scale anthropogenic jet flows (e.g. Matoza et al. 2009;Fee et al. 2010a;Matoza et al. 2013;McKee et al. 2017;Gestrich et al. 2021). Anthropogenic jet noise refers to sounds generated by turbulent exhaust (i.e. a jet flow) exiting aircraft or rocket engines and have been studied in detail for engineering investigations (e.g. Tam et al. 1996;Tam 1998;Tam 2019). Two distinct components of jet noise have been recognised: large-scale turbulence (LST) and fine-scale turbulence (FST; Tam 2019). FST is associated with fine-scale eddies generating acoustics, whereas LST is generated by instability waves forming at the margin of the jet flow (Tam 2019). LST sound radiation is highly directional but may be reduced for volcanoes due to higher temperature jets and a strong diffraction at lower frequencies (Matoza et al. 2009). Jet noise frequency spectra exhibit self-similarity, with spectral shape remaining relatively constant and scaling with frequency, diameter, and velocity (Tam 1998). Studies have shown that this self-similarity can extend to volcanic length scales (metres to hundreds of metres); therefore, similar relationships may exist between volcanic jet noise spectra and vent diameter, jet velocity, and temperature (Matoza et al. 2009;Matoza et al. 2013;McKee et al. 2017;Gestrich et al. 2021). While the spectral shapes are similar, the peak frequencies at volcanoes are generally lower and can be explained through the Strouhal number, a non-dimensional parameter used to describe oscillating flow (Seiner 1984). Quantitative comparisons of similarity spectra were carried out on acoustics recorded during the 2018 Kīlauea fissure eruption which included 80 m high lava fountains (Gestrich et al. 2021). Changes in misfits between the recorded acoustic spectra versus the similarity spectra coincided with variations in eruption dynamics, supporting possible quantitative estimation of eruption flow features McKee et al. 2017;Gestrich et al. 2021).
In this paper, we describe acoustics recorded during the 2021 fissure eruption of Fagradalsfjall, Iceland. In particular, we focus on repetitive lava fountaining that characterised the eruptive activity in May. Each lava fountain was defined by a sequence of activity that can be described by jetting followed by distinct Strombolian explosions (see movie S1 for an example event). Isolation and careful analysis of acoustics during the lava fountaining permits quantification of vent dimensions and highlights the value of acoustic monitoring for future basaltic fissure eruptions.

Fagradalsfjall eruption
The Fagradalsfjall volcanic system is located on the Reykjanes peninsula in south-west Iceland, a 60-km-long transtensional plate boundary between the North American and Eurasian plates (Saemundsson et al. 2020). After a protracted period of earthquake activity along the peninsula since late 2019, an eruption began within the Fagradalsfjall volcanic system on 19 March 2021 (Barsotti et al. 2022;Flóvenz et al. 2022); the first such eruption on the peninsula for nearly 800 years (Saemundsson et al. 2020). The eruption continued until 18 September 2021, at which point the lava flow field had covered an area of 4.8 km 2 with an estimated bulk volume of 0.15 km 3 (Pedersen et al. 2022). Eruptive activity in March and April was spread over 9 distinct vents across a NNE-SSW trending line of 800 m length (Fig. 1a). From late April until the end of the eruption, activity was focused from a single vent (Vent 5 in Barsotti et al. 2022), which we will refer to as the 'main' vent for the remainder of this article. From approximately 01:00 UTC on 2 May, the eruptive activity was characterised by a sequence of repetitive lava fountaining of up to 200 m height above the vent (Barsotti et al. 2022, Fig. 1b;). This period coincided with peak discharge rate of lava from the eruption of 12-13 m 3 ⋅s -1 (Pedersen et al. 2022). Variations in duration, repose intervals, and heights of the lava fountaining were observed throughout May and into early June, eventually transitioning into Page 3 of 18 96 discontinuous lava outpouring that defined the remainder of the eruption (Barsotti et al. 2022). Seismic data analysis (Eibl et al. 2022) and gas emission measurements (Scott et al. in prep.) have suggested the presence of a shallow cavity (< 500 m) beneath the vent where repeated foam collapse in the magma may be driving the distinctive lava fountaining activity.

Acoustic array and processing
The infrasound data presented in this study are from a fourelement campaign array installed on 21 April approximately 800 m north-west of what would eventually be the main vent (Fig. 1a). The 70 m aperture array was equipped with four InfraBSU V2 infrasound sensors, which feature a low noise (5.47 mPa rms at 0.1-20 Hz), flat response from 0.1 to > 40 Hz and inband sensitivity of 45.13 ± 0.23 μV ⋅Pa -1 (Marcillo et al. 2012), along with one co-located 4.5 Hz geophone (not used here), connected to DiGOS DATA-CUBE 3 digitisers (10 V∕

√
Hz analogue-to-digital noise level) recording data at 200 Hz. We estimated a theoretical uncertainty in back-azimuth of 4-5 ∘ and trace velocity of 18-22 m⋅s -1 for this array configuration (Szuberla and Olson 2004). Sensor locations were estimated using time-averaged locations from a handheld GPS receiver. The array ran almost continuously until 30 August, with the exception of 8-17 May when the array had to be removed due to brush fires triggered by the eruption. Technical issues with the sensors from 21 to 23 May led to incorrect values of acoustic pressure levels being recorded, therefore data recorded during this period could only be used for the detection of coherent signals and not for quantifying source parameters. We compare acoustics with seismicity recorded at station MER located north-east of the main vent (Fig. 1a), which used a GEOtiny 10s seismometer recording data at 100 Hz.
Before processing, the acoustic data were filtered with a Butterworth bandpass filter at 0.5-20 Hz. The lower boundary was chosen to remove a strong micro-barom signal that likely obscures signals of interest, and the upper boundary was chosen to ensure as much of the eruptive signal was recorded while removing the majority of anthropogenic noise interference (e.g. helicopters; Barsotti et al. 2022). The data was processed with a least-squares beamforming algorithm (Szuberla and Olson 2004) to identify coherent signals in 10 s windows with 50% overlap. Signals associated with eruptive activity were identified using back-azimuths from 122 to 162 ∘ , trace velocities of 250-400 m⋅s -1 , and a median cross-correlation maximum of > 0.5. To track the number of lava fountaining events occurring during our  Pedersen et al. 2022). Also marked are locations of other vents that were no longer active by May (grey triangles; note that two of the original nine vents merged with or were buried by the main vent), the MER seismic station (yellow square), and the webcam used for estimating lava fountain heights (purple diamond). DEM is from ArcticDEM (Porter et al. 2018 time periods of analysis, we defined an 'event' as a time period when a coherent eruption signal was detected for continuous periods of > 60 s and followed by repose intervals of > 60 s (see Fig. S1 in supplementary material for an illustrative example). For each lava fountain 'event', relative infrasound energy (E a ) was calculated by integrating squared pressure over their whole duration Lyons et al. 2021): We follow the reasoning of Lyons et al. (2021) by using relative acoustic energy because jet noise was a significant component of the wavefield during peak lava fountaining (see below). Therefore, we cannot assume simple spherical acoustic radiation pattern and any measure of energy or intensity on a single array will be relative to itself. Estimating absolute acoustic energy would require observations from multiple arrays at different distances, azimuths, and heights relative to the source ).

Lava fountain height estimations
Lava fountain heights can be estimated from video camera image analysis by developing a kymograph (e.g. Delle Donne and Ripepe 2012; Witt and Walter 2017;Taddeucci et al. 2021). Kymographs were developed for lava fountain events on 5 and 18 May using video recorded by a camera installed on top of Langihryggur approximately 2 km south-east of the main vent by Ríkisútvarpið (RÚV), the Icelandic National Broadcasting Service (Fig. 1a). The camera view angle allows estimations of lava fountain heights below the cone rim due to a cleft in the south-east side of the structure (Fig. 1b, c). Images were recorded at 25 frames per second and lava fountain heights were estimated in each frame relative to sea level. See supplementary information for a full description of how lava fountain heights were estimated as well as a short example video clip from the camera (movie S2). As the camera view angle cannot see the base of the vent (Fig. 1b), we plot all lava fountain heights relative to the lowest section of the breach in the SE section of the cone on 5 May, which was estimated as 238 m a.s.l.

Jet noise similarity spectra fitting
To quantify the presence of jet noise in acoustics generated during lava fountaining activity at Fagradalsfjall, we used an approach recently developed to analyse acoustics during the 2018 Kīlauea eruption (Gestrich et al. 2021). (1) Here we provide a brief description of this approach; see Gestrich et al. (2021) for a detailed method description and discussion. The similarity spectra (i.e. the model spectra) for LST (S LST ) and FST (S FST ) are defined as the following Tam et al. (1996): where f is the frequency, D j is the fully expanded diameter of the jet, r is the source-receiver distance (800 m), and A and B are the amplitudes of the large-scale and fine-scale structures, respectively. The first half the right-hand side of each equation is redefined to C LST and C FST , respectively, so that only one numerical value defines the amplitude of each similarity spectrum. F and G are the spectrum functions of the large-scale and fine-scale structures, respectively, and are dependent on the peak frequency for LST (f L ) and FST (f F ), which ultimately define the characteristic shape of each similarity spectra. In general, LST has a narrower frequency spectrum than FST (Tam 2019).
To fit the similarity spectra with the frequency spectrum of the recorded acoustic data, we use a least-squares method called Trust Region Reflective Algorithm. We use the root mean squared deviation (RMSD) to quantify the difference between the similarity spectrum S and the data spectrum d, defined as Gestrich et al. (2021): for the i th point of a total of n in each spectra. Results are presented within the decibel scale (dB, relative to (20 μ Pa) 2 / Hz). Due to the complexity of the acoustic spectra during each lava fountain event, we focus on using multiple overlapping frequency bands to fit the models to the data and calculate the corresponding root mean squared deviation values. This results in a misfit spectrum for both FST and LST, and subtracting the difference between the two produces a misfit difference spectrogram which shows which of each model is more appropriate for the data over time and different frequency bands (see Fig. S3 in supplementary material for an example). For the jet noise analysis, we limit the frequency bands to the 0.15-20 Hz range using 15 s windows with 90% (2) overlap and we mask out times and frequencies that have an RMSD above 3.5 dB.

Bubble burst modelling
Examinations of the acoustics recorded during the lava fountaining at Fagradalsfjall identified both bubble oscillation ) and Helmholtz resonance (Vergniolle and Caplan-Auerbach 2004) events. Here we describe the method used to describe and quantify the waveforms. A full description of the models can be found in Appendix A.
Waveform inversions were carried out on bubble bursts manually picked from the acoustic record at Fagradalsfjall and we adopted an approach similar to that described for modelling acoustics recorded during the Bogoslof eruption . For bubble oscillation events, we used a model originally formulated by , calculating 10 5 synthetic waveforms from unique combinations of bubble overpressure ΔP (5-80 kPa), initial bubble radius R o (0.1-10 m), and bubble length L (5-200 m), and solved Eq. (A1) with a Runge-Kutta ordinary differential equation solver. For Helmholtz resonance events, we used a model described by , calculating 10 5 synthetic waveforms using Eq. (A6) and directly cross-correlated each with the observed acoustic waveform. The ranges used for the variables were: 50-1x10 4 Pa for ΔP, 0.5-15 m for R hole , and 5-200 m for L. The ranges of each value for each model were chosen based on preliminary testing of the models against selected waveforms. Each synthetic waveform was then directly compared with the observed acoustic waveform and was considered a match if they fell within ± 5% of the observed peak-to-peak pressure and the cross-correlation coefficient maximum was ≥ 0.8. This approach differed slightly for Helmholtz resonance events, as the first (i.e. maximum) peak in the synthetic waveform is an overestimation due to the assumption that the hole radius in the bubble membrane is instantaneously reached . Therefore, we measured synthetic peak-to-peak pressures by ignoring the first 'peak' in the waveform. For each picked bubble event, this method produces a subset of synthetic waveform matches from which we use the median and median absolute deviation values to estimate best-fit bubble parameters and their variability, respectively. This approach results in non-unique solutions for each picked bubble event but these estimations fully explore the model space and may provide insights into changes in bubble source parameters over time.

Results
Due to unforeseen and technical issues (see above), we cannot present a continuous acoustic record of the Fagradalsfjall lava fountain activity through all of May 2021. Here we focus on two time periods and present key observations and analysis from each: 2-8 and 17-24 May; hereafter we will refer to each period as 'early' and 'late' May, respectively.

General observations
Infrasound data recorded during early and late May indicate distinct differences in lava fountaining activity ( Fig. 2) with longer events in the former time period relative to the latter. In general, eruptive acoustics appear with high amplitudes relative to the background activity with peak amplitudes generally occurring at the midpoint or end of each lava fountain event. Acoustic amplitudes in early May generally peak at ≤ 10 Pa (Fig. 2a), whereas amplitudes in late May reached up to and over 20 Pa (Fig. 2c). The frequency content of each fountaining event was generally broadband and may be described as acoustic tremor (Fig. 2b, d) with peak frequencies within 0.5-2.5 Hz.
Array processing results indicate broad fluctuations in lava fountaining activity (Fig. 3). In all, we detected 609 and 1087 lava fountaining events in early and late May, respectively, corresponding to an average lava fountain rate of ∼4.2 and ∼ 6 events per hour. Event durations in early May generally followed a downward trend, falling from < 1000 s down to ∼300 s (Fig. 3a). A break in this trend occurred on early 5 May, when several events occurred with durations up to < 2000 s. Over the same time period, repose intervals followed a generally increasing trend, rising from < 100 to ∼300 s (Fig. 3b). One notable feature is an apparent diurnal variation in repose intervals (and in durations to a lesser degree), with lower values appearing daily at 00:00-06:00 UTC. We observe an apparent correlation with wind velocities recorded at the eruption site ( Fig. S2 in supplementary materials) suggesting reduced signal-to-noise ratios due to wind noise; lower signal-to-noise ratios increase the apparent time interval between lava fountain events. For late May, lava fountaining event durations are shorter and are somewhat more stable with only a slight downward trend from 200 to 100 s (Fig. 3d). Repose intervals between events are longer, fluctuating between 300 to 400 s length, but show no clear indication of any diurnal fluctuations (Fig. 3e). Larger acoustic amplitudes in late May (Fig. 2c) have likely reduced the impact of wind noise on signal-to-noise ratios. Estimated relative acoustic energies during early May show a rising trend in energies over time, from an average of 6.7 Pa 2 ⋅s to 25.8 Pa 2 ⋅s (Fig. 3c). This latter value was sustained on late 17 May, but power levels became more erratic before technical issues affected the sensors at the array (Fig. 3f).
Higher relative acoustics energies of up to 50 Pa 2 ⋅s were estimated after the sensor issues were resolved.
Close inspection of acoustics recorded during individual lava fountaining events and direct comparisons with seismicity generated by the same events indicates a complex activity sequence defined by the superposition of tremor and discrete transients (Fig. 4). In both early and late May, acoustics and seismicity generated during lava fountain events appear almost simultaneously (Fig. 4a, c, e and g). Due to the emergent nature of their onsets, it is challenging to quantify whether there is any significant time lag between the acoustics and seismicity. The first stage of each fountain event is defined by a sustained acoustic tremor with frequencies peaking at < 2.5 Hz (Fig. 4a, b, e and f). As noted previously, the highest acoustic amplitudes are associated with impulsive events with peak frequencies ranging up to 3 Hz; we term these events 'bubble bursts' and describe them in detail below (see below). One distinct difference between events during early and late May is the transition from the tremor to the bubble bursting. Acoustics in early May indicate a more sustained period of tremor for several minutes with multiple concurrent high-amplitude bubble bursting events, followed by a rapid decrease in tremor and a sequence of impulsive events with generally decreasing amplitudes (Fig. 4a, b). In late May, the tremor stage is curtailed and is rapidly replaced by diminishing bubble bursting events (Fig. 4e, f). Seismic amplitudes and frequencies display a distinct difference with the acoustics in both early and late May (Fig. 4c, d, g and h). Amplitudes generally reached a peak in the first stage of the fountaining with slightly reduced amplitudes sustained for most of the fountaining event, before concluding with a second peak followed by gradual decrease to background levels ( Fig. 4c,  g). Notably, the second amplitude peak in seismicity appears after the conclusion of the acoustic tremor during each fountaining event. The frequency content of the seismicity during each fountaining event indicates faint harmonic tremor, with peaks appearing at ∼2.5, ∼ 5, and ∼7.5 Hz (Red triangles in Fig. 4d); these bands are not clear during fountaining events in late May. The distinct acoustic 'bubble-burst' events do not appear clearly in the seismic data, although some peaks in the seismic 'coda' in early May could be related to the ground coupling of atmospheric waves. The lack of detectable seismic energy during bubble burst events has been observed at other volcanoes (e.g. Bouche et al. 2010;Kremers et al. 2013).
Further insights into the lava fountaining activity can be gleaned by comparing the acoustic and seismic recordings with the lava fountain heights estimated from webcam images (Fig. 5). To aid the comparison, we have plotted smoothed absolute acoustic and seismic amplitudes using a 10 s rolling window. Lava fountain heights during each time period of analysis generally display two different behaviours. In early May, lava fountaining quickly rises to reach a peak at 100-150 m in the first stage before a period of sustained low level fountaining fluctuating at 20-50 m (Fig. 5a). In late Fig. 4 Detailed view of acoustic (a, e) and seismic data (c, g) and their respective frequency spectrograms (b, f, d, h) recorded during lava fountain events on 5 (a-d) and 18 May (e-h). Acoustic data was recorded by one microphone of the array, and seismic data was recorded by a broadband sensor at MER (Fig. 1). Red triangles in panel d indicate frequency peaks at 2.5, 5 and 7.5 Hz associated with harmonic tremor during this event May, lava fountain heights quickly reach their peak values at 100-130 m but are instead followed by a rapid decrease down to low levels without any period of sustained fountaining (Fig. 5b). In both time periods, fountaining appears to commence almost simultaneously with the emergence of acoustic and seismic activity. However, the peak lava fountaining does not correlate with the highest amplitudes of the seismo-acoustic wavefield. Interestingly, decreases in lava fountain heights at the end of each event are simultaneous with decreases in acoustic amplitudes, but seismicity generally reaches a peak amplitude following a time delay of < 60 s.

Turbulent acoustic tremor
The onset of each lava fountain event in May is defined by the emergence of acoustic tremor with peak frequencies < 5 Hz (Fig. 4a, b, e and f). Here we present results from systematically quantifying the differences between the acoustic spectra and jet noise similarity spectra (Figs. 6 and 7; see Figs. S4, S5, S6 and S7 in supplementary material for more examples). The misfit difference spectrograms indicate relatively strong LST turbulence before and after each lava fountain event (Figs. 6c and 7c). This is to be expected when background noise is typically comprised of turbulence generated by wind (e.g. Raspet and Webster 2015). During each lava fountaining event, we find a generally bimodal distribution in the misfit spectra: FST is more dominant at > 2 Hz, and LST is more dominant below. One notable departure from this distribution is the appearance of more dominant LST spectra at frequencies > 5 Hz centred on the time of peak lava fountain heights during each event (Figs. 6c and 7c). FST spectra remained dominant in a small 'notch' at frequencies from 2-5 Hz.

Bubble burst modelling
The latter period of each lava fountain event is defined by sequences of transient waveforms and interpreted as acoustics from discrete bubble bursts (Fig. 4). Bubble bursts during each time period of analysis were dominated by m-shaped waveforms linked to bubble oscillations (Fig. 8).
In comparison, events with Helmholtz resonance occur relatively rarely but a notable sequence was observed at 05:00 UTC on 5 May (Fig. 9).
Through inspection of the acoustic recordings, we manually picked 105 and 72 bubble oscillation events on 5 and 18 May, respectively. For events picked in early May, peak-topeak amplitudes ranged from 0.42-23.0 Pa with 1.10-4.29 Hz peak frequencies. In comparison, events in late May had peak-to-peak amplitudes ranging from 0.69-36.54 Pa and 0.75-2.13 Hz peak frequencies. Across both time periods, we matched a subset of synthetic waveforms to all but 4 events in late May, with a minimum of 41, maximum of 4402, and a median of 1577 matches (Fig. 8i- Overall, the matched model parameters tended to span the boundary conditions without clustering near end-member values, indicating an acceptable range of initial conditions. However, we find poorly constrained bubble lengths, with median values falling at the midpoint of the boundary conditions with large median absolute deviations (Fig. 8b, f). Nevertheless, we find narrow distributions of initial radii (R o ; Fig. 8a, e) and overpressures (ΔP; Fig. 8d, h), therefore we interpret these parameters as well constrained. In early May, the initial radii range from 0.3 to 3.3 m (Fig. 8a), and overpressures range from 12 to 63 kPa (Fig. 8d). Despite the poor constraint on bubble lengths, we estimate bubble volumes using Eq. (A2) and find a range of 40-3000 m 3    Fig. 8c). By comparison for late May, the initial radii tend to be larger with values from 0.7 to 4.5 m (Fig. 8e), overpressures are approximately similar with estimations from 13 to 57 kPa (Fig. 8h), and volumes are larger at 162-5000 m 3 (Fig. 8g).
As noted previously, bubble burst events displaying Helmholtz resonance waveforms are relatively rare. However, one sequence of bubble bursts at 05:48 UTC on 05 May was notable for their apparent resonant nature with upward frequency 'gliding' (Fig. 9a, b). This particular sequence occurred immediately after a time interval which featured one of the longest fire fountain events and one of the tallest lava fountain events observed (Fig. S10 in supplementary material). We manually picked 14 bubble burst events from the sequence which had 0.48-10.50 Pa peak-to-peak amplitudes and 1.92 to 4.49 Hz peak frequencies. We matched a subset of synthetics to all events, finding 20-1996 matches with a median of 147 (Fig. 9g, h; see Fig. S11 in supplementary material for all events). Similar to results of modelling bubble oscillation events, the matched model parameters tended to span the boundary conditions without clustering near end-member values, indicating an acceptable range of initial conditions. Manual inspection of the matched synthetics relative to the real waveforms finds a good match for the second to fourth peaks after the first, but poor matching thereafter (e.g. Fig. 9g, h). We find narrow distributions of hole radii (R hole ; Fig. 9c) and overpressures (ΔP; Fig. 9f), therefore we interpret these parameters as relatively well constrained. Bubble lengths (L) do not tend to fall at the midpoint of the boundary conditions, but the median absolute deviations are still large (Fig. 9d). We find radii ranging from 1.7 to 4.1 m and bubble lengths ranging from 42 to 163 m. Using Eq. (A19), we estimate volumes of  (Fig. 9e). Lastly, overpressures ranged from a maximum of 3800 Pa down to 150 Pa at the end of the bubble sequence.

Discussion
The lava fountains observed at Fagradalsfjall were exceptional in terms of event lengths and their recurrence intervals. Previously documented lava fountains at Etna and Kīlauea were more prolonged (e.g. up to multiple months; Ulivieri et al. 2013;Ripepe et al. 2018;Patrick et al. 2019;Lyons et al. 2021) and/or had longer repose intervals (e.g. days to weeks ;Parfitt 2004;Vergniolle and Ripepe 2008). Acoustic amplitudes and lava fountaining heights at Fagradalsfjall are on the same order of magnitude as that recorded during the 2018 Kīlauea eruption (Lyons et al. 2021) but lower than those observed at Etna volcano (Vergniolle and Ripepe 2008;Ulivieri et al. 2013). We find that acoustic data recorded during the lava fountains at Fagradalsfjall required multiple analytical tools to help discern the complex sequence of activity during each event. Results from least-squares array processing found nearly 1700 individual lava fountaining events across two time periods of analysis in early and late May (Fig. 3). This amount approximately correlates with analysis carried out using seismic data (Eibl et al. 2022). The duration and repose intervals along with their broad fluctuations over time also correlate with the seismic analysis but we see an apparent diurnal cycle in early May due to lower signal-to-noise ratios from wind noise. Estimated relative acoustic energies in early May show an upward trend (Fig. 3c) that may reflect changes in eruptive activity over time. Analysis of acoustics from the 2018 eruption at Kīlauea has suggested that increases in acoustic energies correlate with increases in lava effusion rates over long time-scales (Patrick et al. 2019;Lyons et al. 2021), akin to what has been previously achieved using seismic data (Eaton et al. 1987;Battaglia et al. 2005). Time-averaged lava discharge rates compiled over the whole eruption suggest there was a rapid increase in late April and early May (Pedersen et al. 2022). Therefore we suggest the increase in relative energy in early May may be an approximate indicator of increases in lava effusion rates at Fagradalsfjall.
Closer inspection of acoustics generated during individual lava fountains suggests a general sequence of activity during each event: acoustic tremor is followed by a transition into and eventually superseded by impulsive bubble bursting (i.e. Strombolian) activity (Fig. 4). The onset of each lava fountain event is not preceded by any apparent seismic or acoustic signals, which contrasts with observations during eruptive activity at Etna volcano (Vergniolle and Ripepe 2008;Ulivieri et al. 2013). Lava fountains from 2007 to 2013 at Etna were preceded by ∼ 90 mins of Strombolian activity with increasingly larger amplitudes (Ulivieri et al. 2013), and no such precursory activity was detected during the lava fountain phase at Fagradalsfjall. It must be noted, however, that the lava fountains at Etna were larger (800 m) and/or longer (up to 1.7 h; Ulivieri et al. 2013) than those observed in Fagradalsfjall. The apparent simultaneous increase of acoustic and seismic activity at the onset of each lava fountain event suggests both wavefields are, at first, generated by activity at or near the surface. However, the wavefields diverge at the end of each lava fountain event when a second peak seismic amplitude appears with a small time lag after acoustic amplitudes decline (Figs. 4 and 5). The lack of higher acoustic amplitudes at this time suggests a subsurface process with no atmospheric coupling. Another possibility is micro-seismicity generated at the margins of the lava field near the vent, similar to that suggested for tremor recorded during the 2014 Holuhraun eruption (Eibl et al. 2017). It is not immediately clear what generated the postlava fountaining tremor process at Fagradalsfjall and further detailed analysis of the seismicity is beyond the ambit of this article. Nevertheless, it is clear that while seismic data can be used to visualise overall trends in lava fountaining (e.g. Eibl et al. 2022), it may overestimate activity durations and underestimate repose intervals due to 'overshooting' of vent activity (Fig. 5). Another notable observation was peak acoustic amplitudes correlating with Strombolian activity instead of peak lava fountaining. Fluctuations in atmospheric pressure recorded by sensors at distance from a volcanic vent have been directly related to acoustic power, which depends strongly on the source radiation mechanism (Woulff and McGetchin 1976). Three mechanisms are recognised: monopoles can be envisioned as an isotropically expanding source (i.e. an explosion), dipoles are equivalent to directional gas flux interacting with solid walls, and quadrupoles correspond to turbulent noise (i.e. a jet engine). Woulff and McGetchin (1976) use empirical constants on the order of 1, 10 -2 and 10 -5 for each source, respectively. Thus, less acoustic power and lower acoustic amplitudes will be generated during the jetting phase (i.e. quadrupole) than the Strombolian phase (i.e. monopoles and/or dipoles) of each lava fountain event.
Multiple recent studies have studied the relationship between acoustic and seismic tremor amplitudes with eruption rates (e.g. Ichihara 2016; Fee et al. 2017;Haney et al. 2018;Sciotto et al. 2019). Pairwise linear relationships between seismic energy, acoustic energy, and magma discharge rates (which can be estimated from plume or lava fountain heights) can suggest a common seismo-acoustic source within the vent or upper conduit (Ichihara 2016;Haney et al. 2018). Changes in the seismic to eruption rate relationship during the course of the same eruption may indicate a change in the source mechanism (Haney et al. 2018). For the six lava fountains plotted in Fig. 5, we find an apparent linear relationship between lava fountain heights with both seismic and acoustic tremors during the first phase of each event (Figs. S12 and S13 in supplementary material). However, a different relationship exists after peak fountaining with diverging behaviour between seismic and acoustic amplitudes relative to the lava fountain heights. As suggested by Haney et al. (2018), this may be the result of a variation in the seismic tremor source mechanism between the jetting and Strombolian phases during each lava fountain event at Fagradalsfjall.
Jet noise similarity spectra fitting highlighted a complex distribution of fine-and large-scale turbulence frequencies during each lava fountain event (Figs. 6 and 7), with the former dominating at > 2 Hz. The complexity of the acoustic wavefield recorded here is further highlighted when compared with the relatively stable and consistent acoustic turbulence recorded during sustained lava fountaining at Kīlauea in 2018 (see Fig. 5 in Gestrich et al. 2021). While the link between changes in spectral properties of jet noise to eruption dynamics is an ongoing area of research (Gestrich et al, Lava Fountain Jet Noise during the 2018 Eruption of Fissure 8 of Kīlauea Volcano, in review), here we see a potential connection with lava fountain heights. We noted the appearance of stronger LST at frequencies > 5 Hz centred on the peak lava fountain heights during each event (Figs. 6c and 7c). LST sound radiation is highly directional, with a narrow coneshaped wavefield pattern centred on the vent, but this cone may be enlargened for higher temperature and higher Mach number jets (Matoza et al. 2009;Viswanathan 2009). The data presented here hints at an interaction between the cone morphology and the resulting turbulent wavefield, but we are wary in presenting an interpretation. Volcanic jetting is extremely complex and poorly understood compared to laboratory jets and modelling the complex wavefield generated during each fountaining event and how it may be affected by the growing cone would require more data, including from dense acoustic sensor deployments as well as high resolution video imagery looking directly into the vent during an eruption.
Inspection of the Strombolian explosions after each lava fountain event in early and late May revealed two kinds of bursting activity: bubble oscillation (Fig. 8) and Helmholtz resonance (Fig. 9). If we assume the R hole parameter for the latter is approximately equivalent to R o for the former, then the values of bubble radii for each type do roughly agree with each other on 5 May. These radii fall within the range of values found for Strombolian explosions at other volcanoes, including Stromboli, Etna, and Shishaldin Vergniolle and Ripepe 2008;Cannata et al. 2009). However, we find a poor constraint on bubble lengths with the bubble oscillation model with median values around the midpoint between bounds used for generating synthetics. If we instead look at the matching bubble length modal values, we find smaller values at 40-70 m which fall into the upper range of values found in previous studies. On the contrary, overpressures are small compared to previous studies, with peak values at least one order of magnitude lower than the minimum values found at Shishaldin . Analytical modelling of Strombolian eruptions finds a nonlinear relationship between gas overpressure within slugs, and the total volume of the slug as well as thickness of magma between the slug and conduit walls (Del Bello et al. 2012). We note that estimated bubble volumes here are also relatively small compared to the previous studies, which may explain the low gas overpressure estimates. Even so, overpressures for Helmholtz resonance events are small relative to those for bubble oscillation events (Fig. 9f). Helmholtz resonance can also occur within an air cavity above the magma (e.g. Fee et al. 2010b;Goto and Johnson 2011) and is not accounted for in the model we use (see Appendix A). Furthermore, the model assumes a cylindrical bubble volume with constant radius which may lead to overestimates in bubble lengths (see Eq. (A19)). Nevertheless, we suggest the relatively large estimated bubble lengths (Fig. 9d) may indicate an unusually deep magma level in the upper conduit and resonance occurred in an air cavity above the Strombolian explosions.
Laboratory experiments have shown that gas slugs occupy a large portion of the upper conduit with a thin magma film between the slug and conduit walls (e.g. James et al. 2006;Del Bello et al. 2012), therefore the maximum estimated bubble radii should be only slightly less than the conduit radius at the burst point (i.e. magma surface). Images of the vent captured by drone on the morning of 5 May suggested an upper conduit radius of approximately 3.5 m (Joël Ruch, pers. communication). Estimated bubble radii for both bubble oscillation and Helmholtz resonance models on 5 May fall beneath or close to this value (Figs. 8a and 9c). We also found larger maximum bubble radii on 18 May, suggesting an increase in conduit radius during the lava fountaining. This is supported by a steady increase in seismic tremor amplitudes during May 2021 which was interpreted as caused by increasing conduit dimensions (Fig. 4c, g;Eibl et al. 2022). Mechanical erosion of upper conduit and vents have been observed to affect acoustic frequencies during eruptions (e.g. Fee et al. 2017;Watson et al. 2020). We observed a decrease in peak frequencies from 1.1-4.3 Hz to 0.75-2.13 Hz between early and late May, therefore we surmise that conduit erosion also occurred at Fagradalsfjall during the lava fountain sequence. In addition, higher acoustic amplitudes were recorded in late May relative to early May (Figs. 2 and 4) and linear acoustic theory suggests a linear relationship between acoustic amplitudes and mass eruption rates (Woulff and McGetchin 1976;Lighthill 1978). Therefore, we speculate that larger amplitudes during shorter lava fountaining events in late May correspond to higher fluxes of gas over shorter time periods (Fig. 2c) due to a wider conduit (Fig. 8a, e).
However, it must be acknowledged that the bubble burst models used here represent first approximations of the source parameters as source-to-receiver path effects due to topography have not been accounted for. Recent work has shown that recorded acoustic waveforms at volcanoes include amplitude distortion as the waves propagate out of the conduit and around topography (e.g. Kim et al. 2015;McKee et al. 2017;Diaz-Moreno et al. 2019;Lacanna and Ripepe 2020). To account for it, 3D finite-difference time domain modelling which requires high resolution digital elevation models should be incorporated into the source parameter estimations (Kim et al. 2015). Digital elevation models were developed for the Fagradalsfjall eruption site on a weekly basis (Pedersen et al. 2022) but changes in the height of the cone between early and late May were only the order of 10s of metres (also evident in webcam images; Figs. 1 and 5). Artificial scaling of digital elevation models generated after the 2018 Kīlauea eruption indicated that changes in cone height during the activity had no significant effect on acoustic waveforms with frequencies < 1.5 Hz (Gestrich et al. in review). Therefore, viable comparisons may be made between parameters derived from acoustic waveforms recorded at Fagradalsfjall in early and late May (Fig. 8). Nevertheless, future acoustic studies of eruptions with rapidly developing cones and/or lava fields should look to incorporate daily or more frequent measurements of the topography between the source and receivers.
Nevertheless, constraints on the upper conduit dimensions from Strombolian explosion modelling can, in turn, be used to help define the Strouhal number St which connects the peak frequency F p , jet diameter D j and jet velocity U j via St = f p D j /U j . Frequencies during peak lava fountaining ranged from 0.5 to 2 Hz, peaking at 1-2 Hz (Figs. 6b and 7b). Lava fountain heights H can be related to gas velocity v at the vent by v = √ 2gH (Wilson 1980). Peak lava fountain heights during early and late May ranged from 110 to 145 m (Fig. 5), giving peak gas velocities (i.e. jet velocities) of 46.4 to 53.3 m⋅s -1 . Together with the maximum bubble radii found in early and late May (Fig. 8), we estimate St of 0.26-0.43 for the lava fountain sequence at Fagradalsfjall. These represent maximum values as we likely underestimate lava fountain heights since the camera image angle could not see the vent base (Fig. 1b) and air friction on ascending lava is neglected. Nevertheless, these estimations fall within the range of St values previously estimated from volcanic jetting (Matoza et al. 2009;McKee et al. 2017), but lower than estimated for Stromboli volcano (1.2-1.8; Taddeucci et al. 2014). However, we note that peak acoustic frequencies at Stromboli were measured up to 305 Hz (Taddeucci et al. 2014), hence the high St values. Finally, we note that peak acoustic frequencies during each event (red dots in Figs. 6b and 7b) weakly correlate with lava fountain heights, suggesting a relatively constant St value on short-term intervals.

Conclusions
Acoustic recordings can provide important constraints on the eruptive dynamics and source parameters of volcanic eruptions, but data from basaltic fissure eruptions with Hawaiianstyle lava fountains are relatively rare. Here we present an account of acoustic data recorded by a sensor array during the lava fountain sequence of the 2021 Fagradalsfjall eruption in Iceland with analysis focused on activity in early and late May. We find that the array successfully tracked variations in event durations and repose intervals over time but was affected by wind noise. Relative acoustic energies of coherent acoustic detections in early May showed an upward trend in values which may be linked to increasing effusion rates. Detailed inspection of acoustics, and comparisons with seismic data and lava fountain height measurements, revealed a dynamic eruptive sequence during each event: acoustic tremor during peak lava fountaining was followed by a transition to Strombolian activity with distinct highamplitude impulsive waveforms. We observe a divergence between acoustic and seismic amplitudes during each event, with the latter appearing to overestimate lava fountain durations. Jet noise similarity spectra fitting finds complex distributions of turbulent acoustics during each event, with changes of spectral shape and frequency correlated with fountain height. Strombolian explosions are mostly composed of events that can be modelled by oscillations of bursting gas slugs at the top of the magma column, with a minor number of events exhibiting Helmholtz resonance behaviour instead. Modelling of each source type finds maximum bubble radii of < 3.5 m in early May and < 4.5 m in late May, suggesting a widening of the upper conduit during the lava fountain activity. Finally, we propose that a wider conduit in late May, together with higher acoustic amplitudes, indicates higher gas fluxes through the conduit which resulted in briefer lava fountain events. Overall, this study demonstrates the value in incorporating acoustic microphones into future sensor deployments at basaltic fissure eruptions in Iceland and elsewhere through providing important constraints on eruption dynamics and upper conduit dimensions that may not be feasible otherwise.

A.1 Bubble oscillation model
The general equation for an oscillating bubble relates the kinetic energy exchange between a thin bubble cap and potential energy of the gas  where R eq and V eq are the bubble equilibrium radius and volume, respectively, h eq is the bubble membrane thickness at equilibrium (0.15 m, based on observations of deposit sizes around the eruption site), P atm is atmospheric pressure (10 5 Pa), ρ l is the liquid density (2700 kg ⋅ m -3 ), μ l is the liquid viscosity (32 Pa⋅s; Soldatti et al. in prep.), and V g is equal to: where R o is the initial radius of the bubble; V g is a function of the dimensionless bubble radius ε. A sensitivity analysis by  suggests values of h eq at < 0.1 m and > 0.3 m produce unrealistic values in overpressure and bubble length, respectively. The value of μ l does not take into account the effects of crystals and bubbles, the presence of which may increase lava viscosity. However, uncertainties (A1) 0 =̈+ 12 l l R 2 eq ̇+ P atm l R eq h eq in liquid viscosity should not play a major role in bubble dynamics at the air-magma surface ). Heat transfer within large bubbles is adiabatic (Plesset and Prosperetti 1977), therefore internal bubble pressure will follow variations in volume with a ratio of specific heat γ, equal to 1.1 for hot gases (Lighthill 1978). Bubble radius R can be described by its variation from the equilibrium radius R eq  where: in which L is the length of the cylindrical bubble within the conduit, and ΔP is the source pressure, or amount the bubble initially at rest at the magma-air interface is suddenly overpressured.
The pressure recorded at the sensor p ac − p air at time t is  where r is the source-receiver distance (0.8 km), v is sound velocity in air (340 m⋅s -1 ), R is the bubble radius, and ρ air is the air density (1.1 kg ⋅ m -1 ). Therefore, calculation of synthetic pressure waveforms primarily depends on only three variables: ΔP, R o , and L. A full technical discussion of the model can be found in .

A.2 Helmholtz resonance modelling
During Helmholtz resonance, sound is produced by the motion of gas rushing from a small tube into an infinite space, where the tube length is equal to the thickness of the liquid layer (i.e. a hole in the bubble membrane). Sound emission is approximated by that of a piston mounted on an infinite baffle, producing sound as a monopole source in the far field (Spiel 1992). For a piston emitting sound in a halfspace : where ξ is the displacement of air. If the resonator dimensions are small compared to the acoustic wavelength, the behaviour of an undriven Helmholtz resonator is: where m helm , κ helm , and s helm are the mass, the resistance coefficient, and the stiffness coefficient of the oscillator, respectively. These are defined as: where S hole is the hole area, V helm is the volume of the resonator, and is the effective length of the orifice. The oscillator effective length is longer than the geometrical length (i.e. the bubble membrane thickness, h eq ) because some of the air beyond the hole also moves (Spiel 1992). The effective length corresponds essentially to end corrections of a flanged piston which are difficult to estimate since the behaviour of gas about the hole is nonlinear . Spiel (1992) suggests a value of 8R hole /3π ≤ ≤ 16R hole /3π, therefore we use = 4R hole /π as a representative value.
A damped harmonic solution can be defined as: where A and ϕ are arbitrary constants to be determined by initial conditions. The relaxation time τ and the radian frequency ω are: where ω 0 is the radian frequency without damping, defined as : The initial force on the mass of resonating air will be ΔPS hole , the initial speed of the mass will be zero, and the mass of the rupturing bubble membrane will be ignored. These conditions lead to: Finally, the following equations given by  are adapted for an assumed cylindrical bubble of length L with a volume V helm : where R hole is the radius of the hole, which is assumed to be set at its value instantaneously and remains constant over time. Ultimately, calculation of synthetic pressure waveforms from Helmholtz resonance primarily depends on only three variables: ΔP, R hole , and L.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s00445-022-01602-3. Figures S1-S13, a detailed description of lava fountain height estimation methodology, and two videos of the lava fountain activity.