Lava fountains are influenced by many different parameters, like viscosity, volatile content, and magma recharge, that affect the velocity of magma and gas before fragmentation, during fragmentation, and when entering the surrounding atmosphere (Parfitt et al., 1995). In low viscosity, basaltic magma typical for eruptions of Kīlauea, exsolved volatiles form bubbles that are buoyant with respect to the surrounding melt and can move vertically through it (Gonnermann and Manga, 2007). If the magma rise speed is relatively slow, the bubbles can coalesce and form large bubbles that burst at the surface of the magma column in distinct explosions, which are classified as Strombolian-style explosion (Parfitt, 2004; Taddeucci et al., 2015). On the other hand, Hawaiian-style lava fountains can occur when the magma rise speed is relatively fast so that the exsolved bubbles cannot coalesce (Wilson and Head, 1981). The transition between Hawaiian- and Strombolian-style eruptions is marked by lower fountain heights and occasional distinct explosions. When magma rises in the upper crust, the decrease in pressure leads to further gas expansion, acceleration and disruption (Parfitt and Wilson, 1995). Close to the surface, re-entrapment of lava can significantly reduce lava fountain height (Wilson et al., 1995). This happens when lava clasts fall back into the vent and get mixed with the freshly erupting magma. The effect gets noticeable when pyroclasts falling close to the vent form a pond surrounding the lava fountain, which can, in extreme cases, suppress a fountain completely (Parfitt et al., 1995). As the ascending porous magma fragments and interacts with the lava pond at the vent, pyroclast ejection velocity is reduced. After entering the atmosphere, without any additional forces, the pyroclasts can be considered to generally obey a simplistic ballistic equation for fountain height H = u ²/(2g), with g being the gravitational acceleration and u the initial velocity of the observed object (Wilson and Head, 1981). These pyroclast can be accelerated by the drag force applied by the gas expansion (Parfitt et al., 1995). Therefore clast and gas velocity above the vent are mutually related (Wilson et al., 1995).

Quantifying lava fountain properties is challenging, particularly in real-time. Video, radar, gas sampling, and other techniques provide insight into fountain dynamics but are hampered by measurement uncertainties and deployment and safety challenges (Evans and Staudacher, 2001; Witt and Walter, 2017; Mereu et al., 2020). Infrasound provides one technique for examining lava fountains, as they can be prolific producers of low-frequency and audible acoustic waves (e.g., Fee et al., 2011; Sciotto et al., 2019; Lamb et al., 2022). Here, we will use the analogy between laboratory jets and the gas release during lava fountain events to determine the gas velocity and vent diameter. This analysis and interpretation has the potential to be used on real-time data in the future, which are available for many volcanoes. The use of real-time infrasound signals, as opposed to visual measurements, has the advantage to work remotely without interruptions due to visibility. We note that the influence of wind and topographic effects on sound propagation should also be considered, and here we address those in the Supplementary Section S1.1 and Sections 4.2, 5.3.

Kīlauea volcano on the Island of Hawai‘i is the youngest and most active Hawaiian volcano with recent eruptions at the summit in the Halema‘uma‘u crater and along the East Rift Zone (ERZ) at Pu‘u‘ō‘ō vent (labeled in Figure 1) and the lower East Rift Zone (LERZ), which is the most eastern part of the ERZ (Neal et al., 2019). After a period of sustained ground inflation, the activity in 2018 accelerated, starting with the collapse of the Pu‘u‘ō‘ō vent on 30 April 2018, accompanied by increased seismicity and the drainage of the Halema‘uma‘u lava lake in the following days (Patrick et al., 2018; Anderson et al., 2019; Neal et al., 2019). On May 3, an eruption in the LERZ began with the opening of the first active fissure (fissure 1) Neal et al. (2019). Eruptive activity was coincident with significant deformation of the volcano, including a M 6.9 earthquake on May 4 and widespread collapse of the summit caldera for the duration of the eruption. After opening multiple fissures in various locations along the LERZ, the activity focused on fissure 8 (F8), building a cone later named Ahu‘ailā‘au, at the end of May. The location and post-eruptive DEM is shown in Figure 1. The eruption rate diminished significantly on August 5, with weak residual effusion until September 5, when the eruption ended.

The activity at fissure 8 primarily consisted of lava fountaining up to 80 m in height (Neal et al., 2019), feeding a spillway that broadened into a large lava flow north and east of the fissure (Dietterich et al., 2021; Lyons et al., 2021). The erupted basalts had similar compositions to eruptions before 2018 that were sourced from the Kīlauea summit region (Gansecki et al., 2019; Pietruszka et al., 2021). The effusion rate grew slowly, peaking in mid-June at approximately 350 m³s⁻¹ and continuing with relatively high effusion rates of ca. 250 m³s⁻¹ until the end of the eruption in the beginning of August (Dietterich et al., 2021). Beginning as early as mid-June, surges in effusion rate occurred, which are described by Patrick et al. (2019) as long-term fluctuations following collapse events at the volcano’s summit. They were characterized by an increase in bulk effusion rate for approximately 4–6 h followed by a slower decrease. A correlation between bulk-effusion rate and infrasound points to a similar gas content and outgassing mechanism throughout (Patrick et al., 2019; Dietterich et al., 2021). Pulses were another short-duration fluctuation in bulk lava effusion rate, presumed to be due to variations in outgassing efficiency (Patrick et al., 2019). They led to a change of the center of the infrasound source between vent location and spillway and occurred in the period of mid July through the beginning of August (Lyons et al., 2021). Due to the large and rapid changes in infrasound location and outgassing mechanisms during this period of pulsing, we will not focus on the time between July 14–31 in our analysis. Gas measurements and estimations show a similar general trend as the bulk effusion rate with an increase in SO₂ throughout the beginning of June with a peak in mid-June (Kern et al., 2020; Vernier et al., 2020; Lerner et al., 2021).

Jetting is the introduction of a fast fluid into a comparably still medium through a small opening. Volcanic eruptions are known to rapidly discharge gas and solids into the atmosphere. This initial momentum-driven phase has been called jetting in the literature (e.g., Kieffer, 1984; Roche and Carazzo, 2019), pointing to the similar fluid dynamics that are observed for anthropogenic and laboratory jets. This analogy is strengthened by the nature of self similarity within jets, which means that similar features are observed at a wide range of jet dimensions, ranging from small laboratory jets with nozzle openings on the order of millimeters to centimeters to volcanic jets that can be up to hundreds of meters (Matoza et al., 2009; Fee et al., 2010b).

The sound produced by jets has been studied extensively by acousticians and engineers with goals ranging from fundamental research to noise reduction and optimized performance of industrial jets (e.g., Seiner and Gilinsky, 1995). Tam et al. (1996) identified two scales of turbulence that produce jet noise with distinct spectral shapes. Large scale turbulence (LST) noise is produced by coherent instability waves that emanate close to the nozzle region. This noise is dominant in a cone-shaped area around the center axis of the jet. In the jet noise literature, observation angles are typically measured with respect to the downward pointing jet axis. This means that 180° corresponds to a position exactly above the jet nozzle. Large scale turbulence is typically observed at angles above 120° (Tam et al., 2008). Fine scale turbulence (FST) noise, on the other hand, is produced by small incoherent turbulent motion and dominates outside of the LST cone. The spectral shapes described by (Tam et al., 1996) have been used to identify low frequency jet noise during volcanic eruptions (e.g., Matoza et al., 2009, 2013; Taddeucci et al., 2014; Gestrich et al., 2021; Lamb et al., 2022). Jet flow can exhibit a variety of other features that produce acoustic signals such as the vortex ring at the beginning of the jet formation (Fernández and Sesterhenn, 2017; Taddeucci et al., 2021) and crackle during supersonic jet flow (Fee et al., 2013). However, due to the sustained nature of the lava fountain activity, we do not expect any vortex ring formation and we will not investigate crackle noise as there is are no audible reports of such and a lack of skewness (see Supplementary Figure S1 and Section 5.2). Other jet specific noises typical for laboratory jets such as broadband shock noise and screech tones have typically been neglected for volcanic infrasound (Matoza et al., 2009, 2013). Laboratory jets used for studies such as Tam et al. (1996) are typically gas-only jets through comparably small nozzles and with a maximum jet temperature T _j of 3.2 times the ambient temperature T _a (Viswanathan, 2004). In contrast, volcanic jets typically comprise a mixture of gas, solid particles and molten rock in various ratios and temperatures of up to 1,150°C (Gansecki et al., 2019), which is 57 times an ambient temperature of T _a =20°C or 2.5–5.75 times an ambient temperature close to a lava fountain of T _s =200–400°C (Namiki et al., 2021). The effect of the solids and liquids and the increased temperature ratios are not well understood yet and will need future investigation.

Matoza et al. (2013) suggested that some scaling relationships developed for anthropogenic jets may be carefully applied to volcanic jets, such as the relationship between acoustic amplitude and jet velocity and diameter. Note that although similarities between volcanic and anthropogenic jets and their associated jet noise exist, notable differences, mentioned above, and inherent uncertainties require careful evaluation of their applicability.

Acoustic amplitude is usually expressed as sound pressure level (SPL), acoustic intensity, I, or acoustic power, Π, which we define following Garcés et al. (2013) and shown in Eq. 1; Supplementary Eqs S1–S5. Definitions and usage of these quantities vary slightly in the literature, especially with regards to differentiating SPL and intensity. In the following we restrict our analysis to relatively robust scaling relationships that have been developed for laboratory-scale jet noise.

The sound pressure level (SPL) is a parameter that represents the acoustic amplitude. Following Garcés et al. (2013) we define the SPL as: S P L = 10 ⁡ log 10 p rms 2 p ref 2 (1) w i t h : p rms 2 = ⟨ p 2 ⟩ = 1 T s ∫ 0 T s p 2 t d t , (2) With p being the recorded pressure, T _s being a time window of choice and p _ref = 2 × 10⁵ Pa the reference pressure.

Previous studies such as Viswanathan (2004), Viswanathan (2006), Viswanathan (2009), and Tam (2019) showed direct proportionalities between the acoustic intensity, power, and the jet velocity, v _j , and diameter, D _j , for anthropogenic jets, with the applicability to volcanic jets summarized by Matoza et al. (2013). The jet velocity dependence has exponents n _I (θ, T _j /T _a ) and n _Π(T _j /T _a ) for the intensity and power and are discussed below in more detail. Viswanathan (2006) and Matoza et al. (2013) expressed the proportionality between acoustic intensity, I, and power, Π, and jet velocity as. I ∝ v j c n I (3) Π ∝ v j c n Π , (4) With c being the speed of sound.

Furthermore, Matoza et al. (2013), following Tam (2005), reported a direct relationship between acoustic intensity (or sound pressure level) and squared jet diameter D _j . I ∝ D j 2 r 2 v j c n I ρ 0 c p 0 2 ︸ atmospheric conditions , (5) with ρ ₀ and p ₀ being the ambient density and pressure, both of which we assume to be constant. The jet diameter is normalized by the squared source-distance r to account for geometrical spreading of the acoustic wave.

Since the SPL is directly related to the acoustic intensity (see Supplementary Eq. S4), we can express the proportionality in Eq. 3 as. 1 0 S P L / 10 ∝ D j 2 r 2 v j c n I (6) ⇔ S P L ∝ 10 ⁡ log 10 D j 2 r 2 + 10 n I ⁡ log 10 v j c . (7) We note that our infrasound data are recorded in the far-field of the source with the distance being approximately r=500 m. The far-field condition is met for r ≫ λ 2 π ≈ 54  m at a frequency of 1 Hz and speed of sound of 335 m/s.

Importantly, we note the jet velocity exponents n _I for the acoustic intensity and n _Π for the acoustic power are different due to the directionality (directivity) of jet noise (Matoza et al., 2013). The velocity exponent for the acoustic power n _Π does not depend on the observation angle because the power is calculated by integrating over all angles. It has a slight temperature ratio dependence and varies between approximately n _Π = 7.98 to 8.74 for temperature ratios T _j /T _a = 3.2 to 1 (Viswanathan, 2004). Overall there is a slight decrease of the exponent value with rising temperature ratio but the exponent generally stays close to 8. The velocity exponent for the acoustic intensity n _I varies strongly with temperature ratio and observation angle. In Figure 16 in Viswanathan (2006) they show that at 90° relative to the jet axis (i.e., horizontal to nozzle) the exponent varies between approximately 5.8–8.1 for temperature ratios T _j /T _a = 3.2 to 1. The closer the angle gets to the jet axis, the more similar are values for n _I for different temperature ratios. The exponents generally increase the closer the angle to the jet axis and get up to approximately 10. Note that observed velocity exponents for jet noise differ from the simplified cases of monopole, dipole, or quadrupole radiation (Woulff and McGetchin, 1976). These cases are not necessarily applicable in the case of sustained jetting as explained in Matoza et al. (2013).

Viswanathan (2009) showed that the jet noise peak Strouhal number stays the same for different temperatures and jet velocities for an observation angle of approximately 90° from the jet axis. The assumption of a constant Strouhal number has been used in numerous volcano studies (Matoza et al., 2010; Fee et al., 2013; McKee et al., 2017; Haney et al., 2018). The non-dimensional peak Strouhal number is most commonly expressed as S t p = f p D j v j , (8) with f _p being the peak frequency (Mathews et al., 2021). Since the Strouhal number is dependent on the velocity and jet diameter, an increase in jet velocity or decrease in jet diameter results in an increase in peak frequency to accommodate the constant Strouhal number.

Uncrewed Aircraft Systems (UAS) videos of the active fissure were taken most days throughout the eruption and published by DeSmither et al. (2021). Each video is only a few minutes long but provides an overall view of the activity during that day. Example videos from June 2, 10, 16, and 24 are provided in the Supplementary Video S1–S4. The videos are mostly oblique flybys, meaning that the angle, distance, and time varies during and between videos. We use these videos for qualitative description of the eruption dynamics and how it evolved through time. Selected single frames showing representative activity for certain days in June are shown in Figure 2 and described in Section 5.1. Ground-based time-lapse imagery was acquired and analyzed to derive outlines of the change of the Ahu‘ailā‘au cone profile shown in Figure 3. The location of the time-lapse camera is shown in Figure 1A.The post-eruptive topography is documented via an airborne lidar survey in July 2019 (Mosbrucker et al., 2020) and the DEM with 1 m resolution is shown in Figure 1A.

We use campaign infrasound data collected during the eruption of Kīlauea volcano in 2018. Four infrasound elements were deployed about 500 m from fissure 8 to form a small (∼ 35 m aperture) array (Figure 1A). Each element consisted of a Chaparral 60 UHP sensor with a flat response between 0.03 and 200 Hz and was sampled at 400 Hz using a DiGOS DATA-CUBE³ digitizer. These data are available at the IRIS DMC under the network code 5L (http://ds.iris.edu/mda/5L/). Previous research using these data includes analysis of changes in acoustic properties during surges and pulses (Patrick et al., 2019; Lyons et al., 2021) as well as an analysis of the jet noise similarity fitting during (June 16) and after (August 5) the eruption (Gestrich et al., 2021).

We process the infrasound data to determine the amplitudes, spectral characteristics, and source location. First, all infrasound data are downsampled to 50 Hz to reduce computation time. We use least squares array processing techniques to calculate the backazimuth of the incident infrasound waves and delay-and-sum beamforming to increase the signal to noise ratio of the waveform (Szuberla and Olson, 2004; Bishop et al., 2020). Waveforms are filtered between 0.2 and 10 Hz before array processing, which is the dominant frequency band of the signal. We compute the backazimuth using 10 s windows with 25% overlap and then calculate the median backazimuth for every hour.

Primary characteristics of sustained waveforms are amplitude and frequency content. The power spectral density (PSD) is the spectral energy distribution per time increment. Here, we use Welch’s method, which calculates the Fourier transform of a waveform for overlapping segments that are then averaged (Welch, 1977). We use 1 min Hann windows with 30-s overlap averaged over 1 hour. The spectrogram shows the PSDs plotted as a function of time.

The sound pressure level (SPL) is calculated following equation 1. The data are filtered in different frequency bands and we calculate one SPL value per frequency band per hour.

During the eruption of fissure 8 the Ahu‘ailā‘au cone was built from spatter accumulation as shown in Figures 1–2. The cone grew mainly in the first half of June due to the relatively high fountain height and stagnated in growth after the fountain height decreased. A sketched outline of the cone profile from one viewing angle is shown in Figure 3 with cone dimensions estimated by ground-based observations. The DEM that was captured after the eruption had ended, in 2019, reveals that the difference between maximum cone rim height and the bottom of the cone is approximately L = 25.5 m, and the diameter is approximately a = 100 m as shown in Figure 4A . The northeast corner of the cone had a clear opening where lava flowed out into a spillway and then toward the coastline.

We considered resonance of the cone as a potential infrasonic source, as this has been shown to be significant at Kīlauea (Fee et al., 2010a; Matoza et al., 2010) and other volcanoes (e.g., Johnson et al., 2018; Watson et al., 2020). The resonance frequency for a cavity with one open side can be estimated using f _R = c/(4L) ≈ 3.35 Hz. The resonance quality factor for volcanic craters can be defined as Q = 4L/(πa ²) + 32L/(3π ² a) (Moloney and Hatten, 2001; Johnson et al., 2018; Watson et al., 2020). For the geometry of the Ahu‘ailā‘au cone we get Q = 0.35, which is very low compared to other geometries such as in Villarica (Q ∼ 1–3 (Johnson et al., 2018)) or Etna (Q ∼ 2–5 Watson et al. (2020)). A low quality factor is indicative of inefficient trapping of resonant waves and strong radiation of energy into the atmosphere. We therefore do not expect any resonance of the cone to be a dominant source of the infrasound we observe.

Topographic features are known to diffract infrasonic waves and change the waveform significantly (e.g., Kim and Lees, 2011; Kim et al., 2015; Lacanna and Ripepe, 2020; Bishop et al., 2022). We conduct a first order analysis of topographic diffraction and propagation effects using the Finite Difference Time Domain code infraFDTD by Kim and Lees (2011). We use a 1.675 m spacial grid spacing and 0.0025 s time steps. Using the guideline of > 10 nodes per wavelength (Hagness et al., 2005), we assume accurate frequency resolution up to 20 Hz. Since the DEM only shows the topography after the eruption, we assume that this is the maximum height of the cone as there are no reported collapse events after the eruption slowed down. Although the cone growth in Figure 3 is not symmetric, we analyze the influence of different cone heights by scaling the topography as shown in Figure 4A . Lava fountains are complex, with degassing and fragmentation happening at various locations and heights (Namiki et al., 2021). To investigate the influence of source location we vary the source height between 0 and 50 m above the bottom of the cone. We note that we only vary one parameter (cone height or source height) at a time and set the other parameter to their default, which is either the maximum topography or source height of 0 m. The source waveform is synthetically generated by taking the jet noise similarity spectra and inverse Fourier transforming the spectra to the time domain using random phases with a uniform distribution. The similarity spectra used for this was originally fitted to the infrasound PSD (Gestrich et al., 2021) and adjusted for the radial amplitude loss. The final source waveform has a duration of 10 s and is shown in blue in Figure 4B .

Fissure 8 produced a form of volcanic jet noise. Gestrich et al. (2021) developed a fitting method to automatically fit the similarity spectra to volcanic eruption data. They identified that the spectrum recorded from the eruption of fissure 8 fits the jet noise similarity spectra well. We perform a similar analysis for the entire fissure 8 time series and refer the reader to Gestrich et al. (2021) for details on the method. We calculate the root mean square deviation, the misfit between the LST and FST similarity spectra and every PSD calculated in Section 4.1.

Here, we introduce a method to estimate the jet velocity and diameter. We derive two main equations that use the acoustic parameters SPL (sound pressure level) and f _p (peak frequency) to invert for the volcanic jet parameters D _j (jet diameter) and v _j (jet velocity). To do this we first formulate the two base equations. 1 0 S P L / 10 = k 1 c , r , … ︸ constant D j 2 v j n (9) f p = k 2 S t p ︸ constant v j D j . (10) These are adapted from Eqs 6, 8 and are expressed as functions that use jet parameters (diameter and velocity) to solve for acoustic parameters SPL (Eq. 9) and peak frequency (Eq. 10).

Since it is impossible to retrieve actual values for some of the volcanic jet paramters, we use ratios to express how the change in the infrasound parameters relate to the change in jet dynamics.

We use Eq. 9 and solve for the constant k ₁ (Eq. 11). Assuming that the constants remain the same for the entire dataset, we know that the right side of Eq. 11 has to be the same for time t ₀ and time t _i as shown in Eq. 12. Reorganizing and renaming variable ratios (Eq. 13) leads to Eq. 14, which relates the change in jet diameter, velocity, and SPL. k 1 = 1 0 S P L / 10 D j − 2 v j − n (11) ⇒ 1 0 S P L 0 / 10 D j , 0 − 2 v j , 0 − n = 1 0 S P L i / 10 D j , i − 2 v j , i − n (12) ⇒ D j , i − 2 D j , 0 − 2 ︸ ≔ D ̂ 2 = 1 0 S P L i / 10 1 0 S P L 0 / 10 ︸ ≔ S ̂ v j , i − n v j , 0 − n ︸ ≔ v ̂ − n (13) ⇒ D ̂ = S ̂ v ̂ − n o r : v ̂ = S ̂ D ̂ − 2 1 / n . (14)

We then apply the same steps explained in the previous paragraph for Eq. 10 to derive Eq. 18, which relates the change in jet diameter, velocity, and peak frequency. k 2 = f p D j v j (15) ⇒ f p , 0 D j , 0 v j , 0 = f p , i D j , i v j , i (16) ⇒ D j , i D j , 0 ︸ ≔ D ̂ = f p , 0 f p , i ︸ ≔ f ̂ − 1 v j , i v j , 0 ︸ ≔ v ̂ (17) ⇒ D ̂ = v ̂ f ̂ − 1 o r : v ̂ = D ̂ f ̂ (18)

In summary, the ratios for the jet diameter, velocity, frequency, and SPL number are now: D ̂ = D j , i D j , 0 (19) v ̂ = v j , i v j , 0 (20) f ̂ = f p , i f p , 0 (21) S ̂ = 1 0 S P L i / 10 1 0 S P L 0 / 10 (22)

Eqs 14, 18 have to be true and solvable with the same parameter set. Therefore we first set Eqs 14, 18 equal with regards to the diameter ratio and solve for the velocity ratio: S ̂ v ̂ − n = v ̂ f ̂ − 1 (23) ⇔ v ̂ = S ̂ f ̂ 1 n + 2 (24) And secondly we set Eqs 14, 18 equal with regards to the velocity ratio and solve for the diameter ratio: S ̂ D ̂ − 2 1 / n = D ̂ f ̂ (25) ⇔ D ̂ = f ̂ S ̂ − 1 / n − n n + 2 (26) Eqs 24, 26 are now equations that relate the change in SPL and peak frequency to a change in jet velocity and jet diameter. Due to our ∼ 90 ° observations angle we assume values for n _I are between 5 and 8 (Matoza et al., 2013).

In Figure 2, we show eight single video frames taken on different days to show representative activity from each video. In the beginning of June the videos show 1 – 3 relatively tall lava fountains ejected from the fissure 8 vent (see June 2–3 in Figure 2 and Supplementary Video S1). These fountains have a center of variably coherent, porous magma and already liberated gas. The pyroclasts are decelerated by gravitational forces until they reverse direction and fall down outside of the central core. This makes it difficult to identify the diameter of the central core. The height of the fountain fluctuates within a couple of meters over very short time scales (2–6 s). In the UAS videos the falling pyroclasts obscure the view of the central core and the dynamics of the outflow or potential ponding. However, it is clear that the spatter of the fountains leads to the growth of the cone surrounding the fountain (Figure 3). The DEM shows that by the end of the eruption the cone measured approximately 25 m in depth and 100 m in diameter. Note that the depth is measured after the eruption, when lava had ponded and solidified in the cone or partially drained out.

UAS videos of June 9–10 show an apparent slight decrease in fountain height and growth in cone height, with the fountain height still being taller than the cone height (see Figure 2 and Supplementary Video S2). The fountain is broader and less distinct and seems to occupy a larger area of the cone center. The profile of the cone growth is shown in Figure 3. After June 16 the videos show that the fountain height is significantly lower, not reaching over the cone top (see Figure 2; Supplementary Videos S3,S4). With a fountain lower than the cone, we assume that the cone did not significantly gain in height after this time. The activity could be described as a roaring lava pond that is highly turbulent, with lava being splashed against the inner sides of the cone. It is clear that the entirety of the inner cone is churning and overturning. Anecdotes describe a “drowning fountain”, which indicates a rise of the ponding and decrease of fountain height. However, it is difficult to see whether the pond height has increased in the videos.

Interestingly, the fountain height, which continuously decreased in height between the beginning of June to the end of June, does not seem to correlate with the discharge rates determined by Dietterich et al. (2021), which peaked on June 16, and Plank et al. (2021), which peaked between June 9–19.

Array processing shows that fissure 8 is the dominant infrasound source during the study period. The median cross correlation maxima (MdCCM) calculations in Figure 1B show very high correlation values throughout the study period, with the backazimuth centered around the median of 121.7° for the time before the pulses start. This backazimuth coincides with the lowest point of the fissure, which is centered in the surrounding cone. The maximum extent of backazimuth fluctuations before the pulses is between approximately 119° and 125°, which is still inside the cone. Since these values correspond to the backazimuth that produces the highest correlation, we interpret these results as that the highest coherent acoustic energy comes from the fissure itself. With the start of the pulses after July 14, the backazimuth results are more scattered and pointing intermittently to lower backazimuth angles, which is consistent with the interpretation by Lyons et al. (2021) that associates the lower backazimuths with increased activity in the spillway during the pulses.

The waveforms, shown as beamformed signal in Figure 5A, are dominated by a continuous signal from the volcano without showing signs of distinct explosions or notable short-term (seconds to minutes) deviations in the amplitudes. The amplitude distribution can be expressed in skewness with a positive skewness meaning that most of the waveform having positive pressure values. Opposed to supersonic jets that show a positive waveform skewness of up to ∼3.5 in Fee et al. (2013), the eruption of fissure 8 does not show significant skewness (values < 0.15), which is consistent with subsonic jetting (see Supplementary Figure S1).

In Figure 5B we show all calculated PSDs as a spectrogram as well as the 10 h rolling mean peak frequency as a purple line. Here, the peak frequency is the frequency of highest PSD amplitude between 0.2 and 10 Hz. We choose this frequency range because frequencies below 0.2 Hz can be highly influenced by wind as shown in the Supplementary Figure S2 and (Fee and Garcés, 2007). Above approximately 3 Hz the spectral amplitude steadily declines toward higher frequencies and we choose an upper bound of 10 Hz. This frequency range encompasses the majority of the acoustic signal and allows characterization of the spectral slope (Matoza et al., 2009; Gestrich et al., 2021). Generally the spectrum is stable over the study period. The peak frequency in Figure 5B is ∼ 0.7–1 Hz before mid-June (June 16; marked with a dashed line and different background colors) and then declines to approximately 0.3 Hz. This change is not as sudden as it seems in the plot, but rather represents the transition of peak frequency from a steady decline of the frequencies around 0.7–1 Hz and increase of frequencies around 0.3 Hz. The frequency drop in Figure 5B represents when the lower frequencies dominate the higher frequencies. This change is also apparent in Figure 6 where the blue line shows the median PSD in the first time window (until mid-June) and the red line shows the median PSD for the second time window (mid-June to mid-July). The earlier PSD shape has its maximum around 1 Hz and the later PSD’s maximum is around 0.3 Hz.

It is apparent in Figure 5A that the amplitudes of the waveforms rises before peaking in mid-June (June 16) and then steadily declining until mid-July. To further quantify the amplitude evolution we calculate the 10 h rolling mean SPL between 0.2 and 10 Hz (black line in Figure 5C). It confirms that the SPL in this frequency band rises by about 4 dB in the first 12 days of recording and then declines by almost 10 dB until mid-July. We also calculate the SPL for two different frequency bands, both according to the peaks before (0.5–2 Hz) and after mid-June (0.2–0.5 Hz) as shown in Figure 6. One can see that the SPL of the upper frequency band (0.5–2 Hz) has a very similar time evolution as the overall SPL (0.2–10 Hz). In contrast, the lower frequency band (0.2–0.5 Hz) is first lower in SPL than the upper frequency band but surpasses exactly at the time that the overall and upper frequency band SPL peak (June 16) before it peaks about 3 days later. This shows that the sudden drop in frequency observed in Figure 5B only corresponds to the time when the lower frequency band has a higher SPL than the upper frequency band, which is a gradual change.

The FDTD modeling results show that there is no significant influence of the cone or source height for frequencies below 1.5 Hz. The waveforms computed for the different model setups at the location of the infrasound array (see Figure 4B) are used to calculate each power spectral density (PSD) shown in Figure 4C . To emphasize the influence of topography and source height we calculate the difference between the PSD of the default model setup (full topography and source at no elevation) and the PSDs with varying setups as shown in Figure 4D . Our modeling only predicts major amplitude deviations for frequencies > 1 Hz with increases in amplitude (green color in Figure 4D) or decreases in amplitude (pink color in Figure 4D). There are some trends such as a decrease of approximately 10 dB at 10 Hz for higher source heights and an increase of amplitude of more than 10 dB above 20 Hz for very low cone heights. The deviations are generally at frequencies higher than the peak frequency calculated in Section 5.2.

The jet noise fitting shows good agreement between the jet noise similarity spectra and the recorded infrasound. The overall misfit is calculated for the frequency band of 0.2–10 Hz to stay consistent with the infrasound processing in Section 4.1. Note that these are different bounds compared to those used in Gestrich et al. (2021). The peak frequency is bound between 0.15 and 8 Hz. The misfit difference spectrogram is calculated using 20 overlapping frequency bands within the bounds 0.2–10 Hz with logarithmic width so that the highest frequency of each individual frequency band is ten times the lowest frequency. Figure 7A shows the result of the overall fitting and we observe very low misfit values, similar to Gestrich et al. (2021). After June 16 the LST misfit is slightly lower. The misfit difference spectrogram in Figure 7B also shows a better fit of the FST similarity spectrum after June 16 for almost the entire frequency band. The better fit with the FST similarity spectrum is due to its slightly wider shape as shown in Figure 6. Also, the lower frequency bound of 0.2 Hz is above the local minima which is at approximately 0.1 Hz and therefore the fitting misses the left side of the similarity spectra. We show in the Supplementary Figure S3 the result when fitting between 0.1 and 10 Hz, which is slightly more scattered but in general very similar to Figure 7.

First, we calculate the SPL ratio S ̂ (Eq. 22) and peak frequency ratio f ̂ (Eq. 21) using the 10 h rolling mean of the SPL and peak frequency shown in Figures 5B,C as time series SPL _i and f _p,i . As a reference SPL, SPL ₀, and peak frequency, f _p,0 we calculate the mean SPL and peak frequency values for the first 2 days of recording before June 7, which are SPL ₀ = 94.29 dB and f _p,0 = 0.86 Hz. The resulting ratios using Eqs. 22, 21 are shown in Figure 8A and are the input parameters for the scaling relationships expressed in Eqs 24, 26. The results are shown in Figure 8B. Before June 16, both the velocity and jet diameter ratio vary but generally stay around v ̂ ∼ 1 and D ̂ ∼ 1. After June 16 the diameter ratio increases up to D ̂ ∼ 2.5 and stays around D ̂ ∼ 2.3 for the rest of June, which suggests that the volcanic jet diameter has more than doubled since the beginning of June. The diameter ratio slightly decreases in the beginning of July to around 1.1 D ̂ ∼ 1.8. The velocity ratio shows an opposite trend as it drops on June 16 to v ̂ ∼ 0.85, which is 85% of the velocity in the beginning of June, and continues to slowly decrease down to v ̂ ∼ 0.7 until mid-July.

The uncertainty in the velocity exponent n = 5 to 8 leads to variations in the velocity ratio of up to 0.1 and for the diameter ratio up to 0.3 but does not change the overall observed trend of both ratios.