Methods of acoustic gas flux inversion—Investigation into the initial amplitude of bubble excitation

In recent years, there has been a growing need to refine methods for quantifying the volume of gas released into the water column from marine sediments. This demand is fueled by a need to better understand the natural carbon cycle and a requirement to measure, monitor, and verify offshore carbon capture and storage (CCS) sites.^1–6 While there are a variety of techniques to quantify gas emissions, including direct physical measurements, many of these are severely limited in terms of applicability in regard to spatial and temporal coverage, as well as cost. Optical techniques are generally accurate for a limited bubble size range (bubbles large enough to be accurately resolved) but limited in the area they can survey and dependent on water visibility.^7–9 Active acoustic techniques, although able to cover a large area, require time-consuming processing to provide accurate quantification data and give only a snapshot measurement.^10–15 Furthermore, the accuracy of the assumptions that underpin the conversion of active acoustic data into gas fluxes is still evolving.¹⁰ Passive acoustic inversion techniques, on the other hand, are low energy and capable of long-term deployments in any level of visibility and also provide information about the size distribution of bubbles released from the seep.^7,16–21

Passive acoustic inversion relies on a physical model of the sound emitted by a bubble as it is formed. Originally, this principle was restricted to measuring the size of individual bubbles in low flux environments, via the so-called signature method, to ensure the acoustic signature of bubbles did not overlap.^22–25 This approach relies on knowing the frequency emitted by a bubble immediately after release, which is directly related to its equilibrium radius.^19,24,26,27

Signal processing methods were developed to allow this technique to extend to intermediate flux rates when bubble signatures begin to overlap.²⁸ For high flux rates, when bubble signatures do overlap to a significant degree, an approach based on the spectral density of the acoustic pressure measured in the far field at a known distance from a gas seep is appropriate.^21,29 Using a model of the spectrum of the signature of each bubble formation event, the measured spectrum can be used to infer the bubble size distribution. However, the energy released by an individual bubble remains poorly characterized,²¹ particularly, how it varies with depth and mode of injection. This quantity directly affects the gas flux estimate, for example, if the model assumes the bubble is a stronger sound source than it is, then the inversion will underestimate the gas flux.²¹ Conversely, if a single bubble emits more energy than predicted by the model, then the passive acoustic inversion will overestimate the gas flux. To date, the use of the spectral method in the field has proven effective, providing continuous estimates of gas flux over extended periods of time validated by intermittent physical and optical measurements. Unfortunately, such estimates have contained large uncertainties.^7,19,30

This study examines the acoustic signature of bubbles released from coarse-grained sediment to improve the spectral inversion approach by examining the energy released by an individual bubble and the relationship between the initial radius of excitation of a bubble and the equilibrium radius of a bubble. It should be noted at this stage, while a bubble released into the water column may not be spherical, the use of a radial displacement as a descriptor is common in literature as a proxy for more complex volumetric descriptions as (a) the zeroth-order spherical harmonic displacement was the only one that changed bubble volume (to first order) and so radiated a spherical pressure wave,²⁷ and (b) because this aligned the formulation with the historical use in the West of representing bubble volume by an equivalent spherical bubble of known radius.³¹ We then use this refined acoustic inversion technique to quantify ebullition from a natural CO₂ seep occurring at a site with coarse-grained sediment and compare the estimate to physical and simultaneous optical measurements.

The following summarizes the mathematical model for the sound from a plume, which forms the basis of the inversion approach. The theory assumes that the acoustic pressure signature from a single bubble is an exponentially decaying sinusoid.²¹ 

For a spherical bubble, immediately after being entrained into the water column, a bubble will undergo damped simple harmonic motion. The radius of a bubble, $R$⁠, over time, $t$⁠, can be described by

where $Re$ denotes the real part of a complex number, $δtot$ is the total damping coefficient, $H$ is the Heaviside step function, $ti$ is the moment when the acoustic signal is first detected, $tϑ$ accounts for the propagation time between perturbations of the bubble and the corresponding pressure signal at the receiver, $ti−tϑ$ is the moment when the bubble begins to oscillate, $R0$ is the equilibrium radius of the bubble, and $Rε0i$ is the initial amplitude of displacement of the bubble wall.²¹ 

The natural angular frequency, $ω0$⁠, of the bubble is given by an extended form of Minnaert's equation,^19,27

where $κ$ is the polytropic index of the gas in the bubble, $p0$ is the ambient pressure outside of the bubble, $pv$ is the gas' vapour pressure, $ρ0$ is the density of surrounding liquid, and $η$ is the liquid's shear viscosity.

The oscillating bubble wall creates an acoustic pressure in the far field at a distance, $r$⁠, given by

The squared magnitude of the Fourier transform of the pressure radiated from a single bubble at a distance, r, is

The spectrum, $S(ω)$⁠, of the measured acoustic pressure from a population of bubbles with size distribution D(R₀), which is defined such that $∫R1R2D(R0)dR0$ represents the number of bubbles generated per second with a radius in the range (⁠$R1,R2$⁠) can be expressed as

To estimate the number of bubbles via the far field acoustic signal, one needs to determine what value of D(R₀) generates the measured power spectrum, $S(ω)$⁠. This can be achieved by discretising Eq. (5) into $Nb$ finite radii bins such that the nth radius bin is centered on the radius, $Rnc$⁠. Assuming that the bins are contiguous, the power spectrum can be given as

where $Ψ(n)$ is the bubble generation rate within a radius bin (i.e., the number of bubbles formed per second within the nth radius bin) such that $Ψn=Rn+1c−RncD(Rnc)$⁠. The spectrum can be evaluated at the $Nf$ discrete frequencies, $ωm$⁠, and organized into a column vector, $S¯$⁠, allowing Eq. (6) to be written in matrix form,

where $Ψ¯$ is a column vector containing the elements $Ψ(n)$⁠, and M is an $Nf×Nb$ matrix with elements $Mm,n=X(ωm,Rnc)2$⁠.

The matrix equation (7) can be solved for the number of bubbles generated per second $Ψ¯$ in a radius bin, provided that $S¯$ and $M$ are known, to yield an estimate of the gas flux. However, while $S¯$ can easily be determined accurately with the use of a calibrated hydrophone, $M$ remains difficult to accurately determine. This is because Eq. (4) contains one important unknown value, the initial amplitude of the displacement of the bubble wall when it begins to oscillate, $Rε0i$⁠. Most work on studying this has assumed a model of the form, $Rε0i=χradR0$⁠, where $χrad$ is a constant, representing the ratio of the initial amplitude of displacement of the bubble radius and equilibrium radius. The earliest estimate of $χrad$ for bubbles of radii close to 0.001 mm measured near a small waterfall suggested that²⁴ $χrad=10−5$⁠, and the same data were reanalyzed using a more sophisticated model to yield³² $χrad=10−4$⁠. And, considering data from smaller bubbles (⁠$R0≈$ 0.000 312 m) generated in breaking waves suggested³³ that $χrad=1.5×10−2$⁠. Generating bubbles via sheared flows allowed the number of bubbles (>1000) and the range of radii (⁠$R0=$ 0.0001–0.001 m) measured to be greatly expanded.³⁴ This work suggested that $χrad$ was not a constant and increased with decreasing bubble size. Further, in this data, $χrad$ was seen to vary by up to 2 orders of magnitude. One approach to estimating $χrad$ is to compute the flux of bubbles in an experiment where the flux rate can be controlled and measured directly; such an approach has been found to yield estimates of $χrad$ between¹⁹ $1.4×10−4$ and $5.6×10−4$⁠. The apparent variable nature of $χrad$ across bubble sizes and release mechanisms potentially introduces considerable uncertainty into derived flux estimates.¹⁹ 

None of these estimates for the initial amplitude are based on bubbles emerging from sediments. Further, these estimated values are based on a small set of environmental conditions, for example, they are the range of bubble sizes and generation mechanisms.¹⁹ Consequently, there is a need for calibrated $Rε0i$ values, particularly in conditions replicating those found in the field. This paper undertakes such an exercise.

The experimental evidence supporting the assumption that $χrad$ is a constant is rather weak. Herein, the general approach based on assuming $Rε0i$ is a function of $R0$⁠, i.e., $Rε0i=f(R0)$⁠, is adopted so that Eq. (4) becomes

In this study, we analyze the acoustic emissions from 368 bubbles released from sediment into the water column to determine the $Rε0i$⁠, measuring the frequency of each bubble oscillation alongside the maximum radius of excitation and total energy released.

To measure and analyze the acoustic signature of bubbles released from sediment, air was injected into the base of ∼0.30 m of artificial sediment with resulting emissions into the water column recorded with a calibrated hydrophone [Reson TC 4032, Teledyne Technology, Alton, with a sensitivity of 167 dB re 1 V/μPa at 1 m; 0.30 m above the sediment surface; Figs. 1(a) and 1(b)]. Synthetic sediment (Ballotini; 0.009 m diameter, pebbly gravel on the Wentworth scale) was used to ensure uniformity between the grains and simplify porosity and permeability calculations [Fig. 1(c)]. The experiment was designed so that ebullition was in the center of a large 8 m × 8 m (× 5 m deep) freshwater tank [Fig. 1(a)] to reduce the interference of wall reflections, allowing accurate characterization of the acoustic signatures [Fig. 1(d)]. The gas injection rate at the base of the sediment was adjusted until gas migration through the sediment caused bubbles to enter the water column at less than 1 bubble per second. Bubble radii were determined directly from optical measurements (Sony FDR3000 camera, Minato, Japan; 0.50 m from ebullition site).

The optical footage was used to identify periods where single, discrete bubbles entered the water column to ensure that the acoustic signature of an individual bubble release was analyzed. The acoustic and optical data for each bubble were jointly analyzed to determine the start and end time of the bubble oscillation, as well as the maximum pressure of the bubble signal (⁠$Pmax$⁠). The initial excitation radius of the bubble (⁠$Rε0i$⁠) is then calculated using the pressure signals via ^21,34

where $Pmax$ is the maximum pressure signal recorded at the hydrophone as a result of the bubble oscillation.

The center frequency of each bubble signature was determined by iteratively fitting an exponential decaying sinusoid to the signal based on least squares error criterion. The frequency of this wave was used to calculate the equilibrium bubble radius, $R0$⁠, using Eq. (2). The total energy of each bubble oscillation was determined by integrating the acoustic pressure squared as recorded at the hydrophone between the start (⁠$t1$⁠) and end (⁠$t2$⁠) of the event. This is corrected for range, assuming spherical spreading.

where $c$ is the speed of sound in water (1500 ms⁻¹) and $fs$ is the sampling frequency (96 kHz). To account for background noise (electrical and acoustic), the ambient noise on the measurement was calculated using Eq. (10) with the average energy measured over 1 s being equal to 6.2 × 10⁻⁹ J. Subsequently, 6.2 × 10⁻⁹(⁠$t2−t1$⁠) J were subtracted from each event so that only the energy of the bubble oscillation was observed.

To validate these observations, acoustic data collected from a natural seep in offshore Panarea (Bottaro Crater), Italy, is inverted and the estimated flux compared to that obtained by physical bottle measurements collected by a diver and estimates made from simultaneous optical measurements. Panarea is a small Aeolian Island in the southern Tyrrhenian Sea, ∼20 km SW of the active volcano Stromboli (Fig. 2). As a result of the underlying magma chamber, the waters surrounding the island are host to numerous natural CO₂ seeps.^35–38 The calm clear shallow waters surrounding Panarea make it an ideal natural laboratory for testing ebullition detection and quantification techniques.^39–41 Bubbles from the target seep were released from coarse-grained sediment at a rate of approximately 5 bubbles per second, directly comparable to our tank experiment. Similarly, bubbles released from the seep ranged in size from 0.001 to 0.008 m comparable to the 0.003–0.012 m observed in the tank [Fig. 2(c)]. The seep was, however, located at a depth of 12 m, compared to the depth of 2.5 m in the experiments.

Analysis focusses on 10 min of simultaneous optical and acoustic data acquired using an open frame lander, incorporating a calibrated hydrophone (Geospectrum M36, GTI, Dartmouth, Nova Scotia, Canada), and two inward facing Sony FDR-X3000 cameras, focussing on a single central point with scale boards positioned directly opposite. Optical data were processed using techniques described in Li et al., where the radius of each bubble was calculated from the assumed volume of revolution.⁷ Acoustic data were analyzed in 30 s blocks while the optics was analyzed in 15 s blocks.

In total, 368 single bubble releases were analyzed, and the results are discussed below. Natural frequencies were detected between 300 and 1000 Hz, corresponding to bubble radii between ∼0.003 and 0.012 m (Fig 3).

It is useful to examine the energy released by a bubble oscillation as the degree of variation gives us an insight into how consistent the process is. The energy emitted during a bubble's oscillation varies across nearly 3 orders of magnitude, $10−14$–$10−11$ J, with a mean value of $7.5×10−13$ J. Plotting the energy released against the equilibrium bubble radius [Fig. 3(a)], we see a large degree of scatter for $R0$ < 0.006 m, decreasing for $R0$ > 0.006 m.

The initial amplitude of excitation of a bubble ranged mainly between $10−7$ and $10−6$ m with a mean value of $6.0×10−7$ m. Plotting $Rε0i$ against the equilibrium bubble radius [Fig. 3(b)], we see a strong clustering of the data for larger bubbles (⁠$R0$ > 0.006 m) but, once again, a larger degree of scatter for smaller bubbles (⁠$R0$ < 0.006 m). Statistically, the total correlation between $R0$ and $Rε0i$⁠, here, has a Pearson's coefficient of 0.26, a p-value of <0.001 with a trend line for $R0$ and $Rε0i$ expressed in m, determined using total least squares regression (least absolute residual) of

which has a coefficient of determination of 0.99. While this relationship is highly effective for bubbles $R0$ > 0.006 m with the 10th and 90th percentiles of the estimates being within −10% and +15% of observed values. Below $R0$ < 0.006 m, there is a greater variability in the measurement points about the regression line with a root mean square error of $2.6×10−7$ m compared to an average value of $6×10−7$ m.

There is a weak but statistically significant, positive correlation between $Rε0i$ and $R0$ in the results here, particularly, at larger bubble sizes. This contrasts with prior work which assumes that $Rε0i=χradR0$ for a constant $χrad$⁠. It is important to note that the relationship developed here is an empirical one, which is not based on a physical model, and should only be applied to bubble sizes outside the range tested with care. To test the validity of this relationship for bubbles smaller than those tested here, the results are plotted alongside those of a previous study^19,32,34 [Fig. 3(c)], where observed bubbles ranged in size from ∼0.0001 to 0.001 m. Note that this is not a like for like comparison as the experiments use different generation mechanisms: sheared bubbles in a hydrodynamic current at 2 m water depth and releasing bubbles from sediment. Despite the methodological differences, the trend line derived passes through the center of the data cloud for the previous experiment and is in general good agreement down to the smallest bubble sizes, suggesting that this model is applicable across a wider range of bubble sizes than just those studied herein.

Previous assumptions about $Rε0i$ were made by Bergès et al.¹⁹ who used the 25th and 75th percentile of the $Rε0i$ observations of Deane and Stokes³⁴ to invert gas flux from needles and porous stones in a test tank. Comparing our estimates of $Rε0i$ to those used by Bergès et al. when inverting gas flux in Fig. 3(b), we note that their assumption is consistent with our $Rε0i$ observation between $R0=$ ∼0.0001 and 0.003 m, where this is below the range of most bubble sizes observed here. This is consistent with the previous results, only suggesting that the error bounds provided were more conservative than necessary. However, as (1) bubbles were released via a different mechanism, meaning that $Rε0i$ can be expected to vary slightly; and (2) the bubble sizes ranged between 0.001 and 0.004 m radius, meaning that the majority of $Rε0i$ calculations were consistent with this study [Fig. 3(c)].

This highlights the need for researchers using acoustic flux inversion techniques in the field to compare their natural site to that used to calibrate. How similar the range of bubble sizes, flux rates, release mechanisms (needles, sediment type, etc.), depths, or even free field conditions are between the two will determine how accurate the $Rε0i$ estimates will be. Here, an empirical linear relationship between $Rε0i$ and $R0$ has been established for these specific experimental conditions (bubbles released from course sediment in fresh water at 2.5 m depth at a slow rate of a few bubbles per second). The limits of the conditions under which this relationship can reasonably predict the bubble excitation have not been established.

While our data shows a greater variability in $Rε0i$ at smaller bubble sizes, this has little impact on the accuracy of the passive acoustic inversion technique as larger bubbles (for whom we are able to more accurately predict $Rε0i$⁠) dominate the total volume of gas emitted (i.e., 1000 bubbles with a radius of 0.0006 m contribute the same volume of gas as one bubble with a radius of 0.006 m). Hence, in the field, it is more important for flux inversion purposes to accurately calculate $Rε0i$ for large bubbles.

To test the effectiveness of the passive acoustic inversion technique based on Eq. (11), the flux from data collected in Panarea (see Sec. II B) was compared to that recorded by other techniques. A preliminary analysis of the acoustic data revealed a strong continuous signal below ∼500 Hz (Fig. 4). As optical footage suggested very few bubbles larger than 0.006 m in radius (Fig 5), $f0$ = 540 Hz were released, we consider this low frequency noise to likely be the result of moving sediment grains and/or distant seeps.²⁰ Hence, the acoustic inversion is performed between 600 and 2500 Hz (corresponding to bubble radii between 0.001 and 0.005 m) to avoid interference from non-bubble related signals. The upper and lower estimates of our refined acoustic inversion were created assuming a 15% error in $Rε0i$ based on the 95% confidence interval.

The acoustic inversion technique employing Eq. (11) produces a bubbles size distribution (modal bubble radius 0.0047 m) similar to that produced used by Bergès et al. (modal bubble radius 0.0044 m)¹⁹ with slightly more small bubbles reported by the old inversion technique (Fig. 5). Optical inversion of the Panarea data measured a single peak in the bubble size distribution at 0.0034 ± 0.0006 m.

The flux results of our refined inversion technique using Eq. (11) are presented in Fig. 6 alongside prior inversion results, assuming that $χRad$ is between $1.4×10−4$ and $5.6×10−4$⁠. Our refined inversion technique estimate of the gas flux is 3.2 ± 0.6 L/min, significantly greater than that predicted by the previous acoustic inversion technique (0.8 ± 0.7 L/min) but consistent with the mean flux value derived via optical inversion (3.0 ± 0.8 L/min). Note that the uncertainty in the method based on Eq. (11) has reduced by 78% over that employed by the previous method (reducing from 88% to 19%). However, the gas flux estimated by the refined acoustic inversion is slightly higher, 5%–8%, than that obtained by the physical measurement taken prior to deploying the lander (2.3 L/min). The failure of the prior inversion technique to match other flux estimates here is due to the larger bubble sizes seen in Panarea, since it is based on extrapolating $R0=$ 0.002 m to larger bubbles [Fig. 3(d)].

The overestimate in the acoustic inversion estimates of flux performed using more accurate values of $Rε0i$ could be the result of, at least, two possible factors: (1) the role of seabed and sea surface reflections and (2) secondary bubbles oscillations after release.^17,24 To address the first possibility, a model of the direct acoustic path to the hydrophone along with the path reflecting off the seafloor to the hydrophone was formed. This suggested that the energy of the reflected signal is <10% of the direct path, meaning its contribution to the overall overestimation of the flux is too small to explain the discrepancy. This would suggest that the most likely explanation for our overestimate is additional bubble oscillation occurring within the water column as bubbles merge or fragment. Wide-area video footage at Panarea shows multiple bubbles fragmenting in the water column, and this evidence is reinforced when consideration is given to the comparison of the bubble size distribution determined by optical and acoustic techniques. It may be that large bubbles are released at the seafloor where they are recorded by the hydrophone but not the optics, as they are out of thr camera's field of view. By the time the bubbles rise ∼0.75 m into the view of the camera, some have fragmented into smaller bubbles, and these smaller bubbles are measured by the optics for the first time and the acoustics for a second time, resulting in the secondary acoustic peak alongside the optical peak at ∼0.0035 m. In summary, further developments are needed in the acoustic flux inversion code to account for this complex bubble behavior.

As the precision of acoustic inversion techniques improve it will become increasingly more important to quantify the effect of secondary oscillations. To do so, researchers will need some gauge of how frequently bubbles from the seep are fragmenting or merging as compared to the rate of release. This almost certainly varies from seep to seep based on the number of factors, such as the bubble size distribution, where larger bubbles are more likely to fragment,⁴² and the rate of release itself, with bubbles released in rapid succession being more likely to merge together.¹⁷ One could attempt to adapt breakup frequency equations for some indication of the rate of secondary oscillation, however, these typically assume turbulent flows, unlike those seen at the seep here, and does not account for bubbles merging.⁴³ Assuming that secondary oscillation is solely responsible for our overestimate of flux, then the resulting scale factor would be between 0.6 and 0.9.

To improve quantification of gas flux from natural seeps and anthropogenic sources by passive acoustic techniques, we have investigated the initial amplitude of the excitation of bubbles, $Rε0i$⁠, released from coarse sediment at 2.5 m water depth, a previously unknown parameter in the inversion process.

We have empirically determined a relationship between $Rε0i$ and the equilibrium radius of a bubble, $R0$ (in m), of the form $Rε0i=9.608R0+0.447×10−6$⁠.

Our data show a weak positive correlation, especially at larger bubble sizes, and appears generally consistent with prior investigations of smaller (<1 mm) bubbles.

We demonstrate the potential of our refined inversion technique by inverting the flux from a natural CO₂ seep in Panarea (Italy) and comparing the results to simultaneous optical measurements. Here, we find that after filtering out low frequency, non-bubble related noise, our acoustic flux estimate is consistent with optical measurements and only 5%–8% larger than prior physical measurements. We identify the post release fragmentation of bubbles as a likely explanation of minor overestimation of flux.

The 78% improvement in the precision of the acoustic inversion technique demonstrated here will allow for the more widespread adoption of passive acoustic techniques for the long-term observations of natural seeps and the measuring, monitoring, and verification (MMV) of secure CCS sites.⁶ 

Funding was provided by the European Union's Horizon 2020 research and innovation programme under the Grant No. 654462 (Strategies for Environmental Monitoring of Marine Carbon Capture and Storage). We would like to thank Michele Deponte, Emiliano Gordini, Diego Cotterle, Valentina Esposito from the Italian National Institute of Oceanography and Applied Geophysics for their logistical support and the scientific divers Andrea Fogliozzi and Martina Gaglioti. We are grateful to the European carbon dioxide capture and storage laboratory infrastructure (ECCSEL)—NatLab Italy of Panarea, an ECCSEL-European Research Infrastructure Consortium infrastructure of Istituto Nazionale di Oceanografia e di Geofisica Sperimentale, without whom our field work would not have been possible. We thank John Davis for his contributions to the lander design and fieldwork.