## Abstract

In recent years, because of the importance of leak detection from carbon capture and storage facilities and the need to monitor methane seeps and undersea gas pipelines, there has been an increased requirement for methods of detecting bubbles released from the seabed into the water column. If undetected and uncorrected, such leaks can generate huge financial and environmental losses. This paper describes a theory by which the passive acoustic signals detected by a hydrophone array can be used to quantify gas leakage, providing a practical (as opposed to research), passive and remote detection system which can monitor over a period of years using simple instrumentation. The sensitivity in detecting and quantifying the flux of gas is shown to exceed by more than two orders of magnitude the sensitivity of the current model-based techniques used commercially for gas leaks from large, long pipelines.

## 1. Introduction

In 2007, Leifer & Tang (2007) published records of the passive acoustic emissions from single bubbles emitted from a natural marine hydrocarbon seep. The sparse bubble generation from this modest seep allowed the emission from each bubble to be separately observed, and hence (through comparison with photographic data) tests could be made of the so-called ‘Minnaert formula’ by which the bubble radius can be inferred from the frequency of its ‘signature’ passive acoustic emission. Larger gas emissions (from seeps, ruptured gas pipelines or leaks from facilities for carbon capture and storage) have the advantage of generating stronger signals, such that they can be detected at greater ranges (and so allow a fixed hydrophone array to monitor a larger area continuously for such events). However, these larger leaks come with the disadvantage that the bubble emissions overlap, confounding attempts to count bubbles by identifying individual bubble ‘signatures’, requiring an alternative that can cope with such overlapping. This paper theoretically tests and quantifies such an alternative, which has the potential to estimate the gas flux and bubble population emitted by such large events, and explores the limitations of the approach.

The topic of gas bubbles in marine sediments is of increasing importance (Judd 2003), first, because of the impact those bubbles have on the structural integrity and load-bearing capabilities of the sediment (Wheeler & Gardiner 1989; Sills *et al.* 1991; Briggs & Richardson 1996); second, because of the impact which the bubbles have on any acoustic systems used to characterize the sediment (Anderson & Hampton 1980*a*,*b*; Karpow *et al*. 1996; Anderson *et al*. 1998; Boyle & Chotiros 1998; Wilkens & Richardson 1998; Gardner 2000; Leighton 2007; Leighton & Robb 2008); and third, because the presence of bubbles can be indicative of a range of biological, chemical or geophysical processes. This includes the assessment of gas reserves for fuel and the climatologically important flux of methane from the seabed to the atmosphere (Judd 2003). The ability to generate ‘greenhouse’ warming per molecule of methane gas is at least 20 times that of each CO_{2} molecule (Khalil & Rasmussen 1995), and any assessment of marine gas reserves should also factor in the potential contribution from hydrate dissociation, which will be promoted through warming associated with climate change (Westbrook *et al.* 2009). The assessment of Dillon (1995) is that the global reserve of methane in the form of hydrate is more than twice the worldwide amount of carbon to be found in all known conventional fossil fuels on Earth.

Furthermore, in recent years, proposals that large-scale geological carbon dioxide sequestration be a major component of the global carbon mitigation strategy have given impetus to finding methods for a long-term leak monitoring. The UK seabed capacity for carbon storage is estimated to equate to 100 years of the current power sector output (Gough & Shackley 2005), with the current Government providing up to £1 billion for the first commercial-scale carbon capture and storage (CCS) project, with a further commitment to provide public sector investment in three additional commercial-scale CCS projects. The industry is looking for systems to detect leaks from these facilities. Such leaks would compromise them, and could potentially affect marine ecosystems, ocean acidification and the carbon cycle. However, a successful leak-monitoring solution could also be applied to natural methane seeps and commercial subsea gas pipelines. The methods currently in use for monitoring subsea gas pipeline leakage are not sufficient for the increasingly stringent requirements of regulatory authorities across the world, in terms of sensitivity, ability to localize the leak, and the power and maintenance requirements (Theakston 2004; McStay *et al.* 2005). Sensitivity is a particular issue, since if it is insufficient a leak that remains small may persist for long periods with accumulating effects (financial loss, environmental damage, etc.), and a leak which is growing will not be stopped at an early stage. Gas leakage is a significant issue in an industry where an offshore development can have a $10 billion value (Devold 2009; Det Norske Veritas 2010). Although the authors could find no mention of a passive acoustic detection system for gas pipeline leakage in the peer-reviewed literature, the trade literature (Barbagelata & Barbagelata 2004) describes a system which identifies leaks from the passive acoustic emissions, but cannot quantify the leak, as an inclusion of the approach outlined in this paper would do.

Because of the requirement for practical (as opposed to research) systems to provide sufficient coverage over space and time for detecting bubbles released from the seabed into the water column, there is a particular need for remote and passive methods which can continuously monitor over periods of years using simple instrumentation and provide a quantitative estimate of the gas flux. This paper describes the signal interpretation needed to achieve this by passive acoustic monitoring using a hydrophone array (with a clear potential for transposition to an underwater or surface vehicle), so providing a technique for detecting and quantifying leaks from natural and man-made underwater structures such as methane seeps, CCS facilities, gas pipelines, riser pipelines and underwater hydrocarbon production templates.

Section 2 considers the ‘forward problem’ (predicting the acoustic emission from a bubble released during a leak) and §3 will consider the ‘inverse problem’ (estimating the flow rate in the leak from the measured acoustic signal).

The injection of a single bubble of radius greater than a few tens of micrometres in water produces a ‘signature’ (an exponentially decaying sinusoid at the so-called Minnaert frequency, which simply relates the frequency of the bubble emission to the bubble size). Identification of this signature had, since 1933, been used to assess the size of individual injected bubbles in tanks by a set of authors (Strasberg 1956), and in the 1980s was first used to size a naturally generated bubble population (Leighton & Walton 1987), through the study of waterfalls and streams. The technique was applied to other natural phenomena, e.g. wavebreaking and rain at sea (Updegraff & Anderson 1991; Leighton *et al.* 1998*a*); and methane seeps (Leifer & Tang 2007). The ‘signature’ method was also applied to laboratory scale tests of a range of natural phenomena, including bubbles generated by rain over water (Pumphrey & Walton 1988) and wave breaking (Medwin & Beaky 1989; Kolaini & Crum 1994). However, such techniques become impractical as the bubble signatures overlap, as is the case when the bubble production rate is high (Leighton *et al.* 1991). Nikolovska & Waldmann (2006) proposed extending the detection of single bubble signatures to the realm of larger gas fluxes by first collecting the gas in an underwater inverted funnel, and then releasing the gas from that funnel as large single bubbles, the size of which could be detected by their acoustic signatures. The gas flux estimation would therefore not be based on detection of the acoustic signatures generated on release of the bubbles from the seep/leak, but on release of the gas from a funnel which has been placed over pre-existing gas leaks. While this worked in a small test tank for a funnel of 36 cm maximum diameter, in practical usage it could only quantify leakages over which a funnel had been placed (and would therefore be unsuitable as an alarm for gas pipe leaks). More fundamentally, if the gas has already been collected, then the advantage of releasing it into the water column to detect its volume, as opposed to direct measurement of the volume of captured and retained gas, is not clear: once leaked gas has been collected, it would seem counterproductive to release it in order to measure its volume. However, conceptually this technique does extend the principle for detecting individual non-overlapping bubbles to higher gas production rates, since the bubbles in question were made larger through bubble coalescence in the funnel. Nikolovska & Waldmann (2006) successfully used it in the laboratory to quantify the number and size of air bubbles released, and hence the gas flux, using a single hydrophone, and certainly the technology for seabed gas collection is well developed (Leifer *et al.* 2006). However, the method of Nikolovska & Waldmann (2006) requires that the nozzle chosen to release bubbles from the funnel does so in a manner which allows acoustic identification of single bubble events, which severely limits the scale-up potential. Although the 1 mm internal radius of the nozzle used in their experiment generates single bubbles for small gas fluxes, bubble nozzle dynamics do not scale out of the O(mm) nozzle diameters they used (Clift *et al.* 1978; Leighton *et al.* 1991). In their study, the funnel was 36 cm radius at its collecting base. Using multiple small funnels of this type is clearly not cost-effective, given the number of hydrophones required; and the alternative of scaling up to a large nozzle to emit single separate bubbles from the area required to monitor a large seep field would not work because such nozzles would not emit single bubbles with individually distinct acoustic signatures (Leighton *et al.* 1991). Even over the limited size scale where the technique worked in the laboratory, it cannot measure the bubble size distribution on release from the seep, which is a required input for models of mass transfer in the evolving bubble plume (Maksimov & Sosedko 2001; Maksimov 2003*a*,*b*; Maksimov & Polovinka 2004).

In practical usage, increasing gas fluxes cause bubble signatures to overlap, and none of the above methods become applicable. In such circumstances the individual bubble signatures, and the bubble size distribution must be inverted from the overall spectrum. While there is a call for such inversions in many industrial processes (Boyd & Varley 2001; Manasseh *et al.* 2001), any such inversion which relies on free-field theory (such as that of Minnaert (1933) and Devin (1959)) must be critically examined for errors. This is because while such inversions can always be made to produce an answer, the accuracy of that answer depends on the correctness of the underlying physics, and when bubbles are confined within a vessel, the use of free-field theory can generate very large errors (Leighton *et al.* 1998*b*, 2002). In the type of inversion discussed in this paper, validation of the underlying physics is more difficult than validation of the inversion process itself, as the following sections will demonstrate.

## 2. The forward problem

Consider for the moment a single gas bubble emitted from the leak of interest. The bubble has a vector position *x*_{0} with respect to the origin in a liquid (assumed to be incompressible) of density *ρ*_{0}. In the limit of small amplitude oscillations, and for *kR*_{0}≪1 (valid assumptions for most ocean gas bubbles pulsating at their natural frequencies; Ainslie & Leighton 2009), this approach can be used to find the time history of the pressure measured at vector position ** x**, where

**−**

*x*

*x*_{0}=

**. Assume for the moment that the bubble wall motion is simple harmonic and undamped, at some angular frequency**

*r**ω*. The instantaneous value of the radius is

*R*(

*t*) and this oscillates about the equilibrium position

*R*

_{0}. The wall displacement

*R*

_{ε}(

*t*) oscillates harmonically at a circular frequency

*ω*with an amplitude

*R*

_{ε0}, where

*R*(

*t*)=

*R*

_{0}+

*R*

_{ε}(

*t*) such that . Then 2.1 the negative sign indicating that, in the quasi-static condition, a positive applied pressure will lead to a compression. The factor

*t*

_{ϑ}accounts for the propagation time

*r*/

*c*

_{0}between perturbations of the bubble and the corresponding pressure signal detected at a distance

*r*(Leighton 1994).

In the linear limit the nonlinear terms can be neglected, allowing substitution of the binomial expansion
2.2
into the expression for the oscillatory pressure in the liquid *P*_{b1}(*t*), at a given range *r* from the single bubble, such that the following reduction takes place in the linear limit (Leighton 1994):
2.3
This equation relates to the undamped low-amplitude pulsations at a frequency *ω*_{0} of a bubble in a liquid. In fact the bubble starts pulsating at a certain time, after which those pulsations are damped. Equation (2.3) can be used to predict the pressure ‘signature’ detected by a sensor at range *r* from a single pulsating bubble, provided that there is a suitable model for the time-dependency of *R*_{ε}(*t*) (§5 describes some practical methods of establishing range of the receiver from the bubble). For bubbles entrained in the process of a gas leakage from, say, a carbon storage facility, the following equation provides a useful approximation of the damping-induced time-dependency of this:
2.4
where *H* is the Heaviside step function, *t*_{i} is the moment when the acoustic signal is first detected at the monitor, and *t*_{i}−*t*_{ϑ} the moment when the bubble begins to oscillate. The exponential decay of the bubble, which is pulsating at its natural frequency *ω*_{0} and initial wall amplitude *R*_{ε0i}, is governed by the total dimensionless damping constant, *δ*_{tot} (Leighton 1994, §§3.4 and 4.4.2). The term *ω*_{0} is the natural angular frequency of the bubble, given by
2.5
This adapts the familiar Minnaert frequency to include the effects of vapour pressure *p*_{v}, surface tension *σ* and shear viscosity *η*. The polytropic index *κ* varies between *γ* (the ratio of the specific heat of the gas at constant pressure to that at constant volume) and unity, depending on whether the gas is behaving adiabatically, isothermally, or in some intermediate manner. Note that the use of a polytropic law only adjusts the way gas pressure changes in response to volume changes to account for reversible heat flow between the gas and its surroundings (Leighton 1994, §§3.4 and 4.4.2).

The model of equation (2.5) is only approximate: for example, most acoustic signatures for entrainment show a small ring-up, as well as decay, period. This can be included in equation (2.5) by replacing *H*(*t*−*t*_{ϑ}−*t*_{i}) by *H*(*t*−*t*_{ϑ}−*t*_{i}) (1−*e*^{−Γt}), where *Γ* is a ring-up factor evaluated by empirical observation or from the developing theoretical base (Longuet-Higgins 1990; Pumphrey & Ffowcs Williams 1990; Leighton 1994, §3.7.3; Clarke & Leighton 2000; Deane & Czerski 2008).

Substitution of (2.5) into (2.3) gives the monopole emission detected in the far field from a single bubble:
2.6
with the Fourier transform:
2.7
the squared magnitude of which is
2.8
If the acoustic emissions of the bubbles are all uncorrelated, then the power spectral density, *S*(*ω*), of the far-field acoustic signature of the bubble cloud, at range *r* can be expressed as
2.9
where *D*(*R*_{0}) is the bubble-emission size distribution as a function of *R*_{0}, defined such that represents the number of bubbles generated per second with a radius in the range (*R*_{1}, *R*_{2}). Manipulation of the integral in (2.9) must take into account the dependence of *ω*_{0} on *R*_{0} as expressed through (2.5). Equations (2.8) and (2.9) can together be used to invert measurements of the power spectrum of the sound field from bubble generation through pipe or sequestration facility leakage, as detected in the far field, to obtain an estimate of the rate at which bubbles in given radius size bins are generated. This will in turn give the flow rate.

There is, however, one important unknown in equation (2.8), the initial amplitude of displacement of the bubble wall at the start of the emission as a proportion of the equilibrium bubble radius (*R*_{ε0i}/*R*_{0}). Use of a constant factor would facilitate the inversion (Loewen & Melville 1991), but may not reflect the real relationship. There have been few estimates of *R*_{ε0i}/*R*_{0}. Of those, the majority followed the example of the first (Leighton & Walton 1987) who, rather than applying the model of (2.6), inverted the measured sound following bubble injection using a model of a rigid pulsating sphere for the acoustic emission (where the sphere wall moves but the pressure inside it does not fluctuate, such that the radiated field is based on wall motion only). Using this rigid model, Leighton & Walton (1987) estimated *R*_{ε0i}/*R*_{0}≈10^{−5} for *R*_{0}≈1 mm, and following the same method Pumphrey & Walton (1988) determined *R*_{ε0i}/*R*_{0}≈7.5×10^{−3} for *R*_{0}≈0.44 mm for bubbles entrained by liquid drop impact, and Medwin & Beaky (1989) estimated that *R*_{ε0i}/*R*_{0}≈1.5×10^{−2} for *R*_{0}≈0.312 mm for bubbles entrained under a laboratory breaking wave. However, Leighton (1994, §3.5.1b(ii)) reanalysed the raw data of Leighton & Walton (1987) using the better model (equation (2.6), where changes in bubble volume perturb the gas pressure and hence contribute to the radiated field) to estimate that *R*_{ε0i}/*R*_{0}≈10^{−4} for those injected bubbles. Deane & Stokes (2008) analysed the acoustics generated when 505 bubbles fragmented into pairs in shearing flow, and found that *R*_{ε0i}/*R*_{0} tended to decrease (with some scatter) with increasing bubble radius, from 0.02 for *R*_{0}≈0.064 mm to 0.0002 for *R*_{0}≈2 mm. At any given bubble radius measured, there was around an order of magnitude of scatter in values of *R*_{ε0i}/*R*_{0}. They conclusively identified the mechanism for bubble excitation for the shear fragmentation they measured (which would probably resemble excitation during injection or underwater gas pipe leakages) as a volume change in the bubble gas generated by the rupture of the necking region between the future fragments, the energy being manifest at the moment of fragmentation in the surface tension energy associated with the geometry of the neck (Deane & Czerski 2008; Deane & Stokes 2008). Importantly, they predicted that the forcing function was proportional to the internal gas pressure in the bubble, suggesting the possibility that for a given bubble radius the quantity *R*_{ε0i}/*R*_{0} may be broadly invariant with depth (within the limits of the observed scatter, over much of the bubble size range). There is no published data as to the extent to which this invariance will be observed, or the extent to which these trends for injected and shear-generated bubbles will transpose across the wide range of applications (natural seeps, hydrate dissociation, leaks from pipe and CCS facilities, etc.) considered however. However, the developing theoretical basis supports the position taken here, and a depth-invariant value of *R*_{ε0i}/*R*_{0}=3.7×10^{−4} across the bubble radius range studied here (about 1–5 mm for figures 2 and 4; and 0.2–6 mm for figure 6) represents a starting position which can be modified by radius and depth dependencies should a future extended radius range warrant that. While the developing theoretical basis should provide analytical predictions of *R*_{ε0i}/*R*_{0} in the longer term, calibration exercises need to be undertaken for test leaks typical of those to be found in the field, to determine the only unknown parameter, *R*_{ε0i}/*R*_{0}, (which may depend on *R*_{0}) from data for the correct environmental conditions (especially depth, bubble radius, surfactant presence and bubble production mechanism). There are few controlled studies in conditions relevant to the release of carbon dioxide or methane gas (e.g. liquid half-space above appropriate seabed and a ‘nozzle’, which may be buried or not) to provide testing for models of the acoustic emission under such circumstances, or to provide observations to adjust the models empirically such that the models match the at-sea conditions. The measured damping for free oscillations of propane (and other) bubbles injected into small tanks tended to be considerably greater than that predicted by theory (Leighton & Walton 1987), a feature confirmed by Walton *et al.* (2005) for methane (and other) bubbles, with decay times and quality factors differing from predictions of Devin (1959) by as much as a factor of 2. However, it would be incautious to extend this directly to at-sea conditions since these tests compared injection in environments having high reverberation (which can considerably affect damping, Leighton *et al.* 1998*b*, 2002) with free-field theory. In a rare at-sea study, Leifer & Tang (2007) noted from correlating optical and acoustical measurements that the actual radii of methane bubbles were as much as approximately 20 per cent lower than that predicted by the Minnaert formula. A long-standing departure from Minnaert's prediction is found in the work of Strasberg (1953), who studied the influence of bubble shape or nearby surfaces (such as the seabed or neighbouring bubbles) on the natural frequency. There is therefore potential for the actual emissions to deviate from the standard expressions. These effects might generate significant errors in some inversions, and negligible effects in others (if, for example, a wide distribution of bubble sizes smears out the discrepancy of some bubbles; or if the inversion is for a parameter such as the void fraction, which has such a large potential dynamical range that factors of 2 or 3 can be insignificant for some practical applications).

It has yet to be ascertained to what extent these preliminary findings are typical of the general population of at-sea seeps or leaks from CCS facilities. Certainly, in the oceans natural factors give potentially wide variations in the chemical constitution of the bubble wall and the details of the bubble shape distortion, release and fragmentation (Leighton *et al.* 1991; Maksimov & Leighton 2001, 2011; Winkel *et al.* 2004; Lee *et al.* 2007). These will provide many potential features by which the forward model of equations (2.5) and (2.6) can be adapted as new data come in. It should be noted with caution that the ‘control’ laboratory conditions against which these field observations and models can be compared for validation can also contain features (in terms of reverberation, water cleanliness, etc.) which can generate artefacts (Leighton *et al.* 2002). Therefore, after explaining the method of inversion (§3) and undertaking inversions of constructed datasets (§4), this paper will conclude by addressing to what extent the levels of variability indicated in the preliminary findings of Walton *et al.* (2005) and Leifer & Tang (2007) can affect the bubble size distributions and gas fluxes determined by remote inversion of passive acoustic emissions.

## 3. The inverse problem

In order to estimate the number of bubbles from measurement of the far field power spectrum, one needs to solve (2.9) for *D*(*R*_{0}) given *S*(*ω*). Equation (3.1) shows one way of doing this, discretizing (2.9) into *N*_{b} finite radii bins, such that the *n*th radius bin is centred on the radius and the bin limits set at *R*_{l,n}<*R*_{0}≤*R*_{u,n}, where *R*_{l,n} and *R*_{u,n} are the lower and upper limits of the bin, respectively. In addition, it is assumed that the bins are contiguous, so that *R*_{u,n−1}=*R*_{l,n} for *n*=2,3,…,*N*_{b}, in which case:
3.1
Here, *Ψ*(*n*) is the bubble generation rate within a radius bin (i.e. the number of bubbles formed per second within the *n*th radius bin). This can be expressed as
3.2
In practice, the spectrum *S*(*ω*) is evaluated at a discrete set of frequencies, *ω*_{k}, with the radius bins selected such that, according to (2.5), corresponds to bubbles with an angular resonance frequency *ω*_{k}, in which case, the number of frequencies is equal to the number of bubble radii and one can write
3.3
Here ** Σ** is referred to as the spectral matrix, defined such that it has elements , and

**S**is a column vector containing the measured spectrum

*S*(

*ω*

_{k}). The object of the exercise is to estimate the vector (

*N*

_{b}×1) of bubble generation rates,

**, the elements of which are**

*Ψ**Ψ*(

*n*) as defined in (3.2). An inherent assumption is that the majority of the acoustic energy is emitted by bubbles upon their generation, rather than their re-excitation, which is generally justifiable (Pumphrey & Ffowcs Williams 1990), although there will be double counting of gas if a counted bubble then fragments.

The matrix ** Σ** will be square (

*N*

_{b}×

*N*

_{b}) as long as the number of frequencies is equal to the number of radius bins. In such a case, in principle, the population can be estimated from the measured spectrum simply using matrix inversion, so that 3.4 The spectral matrix is not symmetric, and furthermore it is frequently ill-conditioned. Inevitably, the measured spectral levels in

**S**will include some noise contributions, which the ill-conditioning of

**can magnify to a degree that leads to unphysical results, i.e. negative estimated bubble generation rates. Thus, it is advisable to apply some form of regularization to the solution. One simple form of this is called the Tikhonov regularization, also known as ridge regression (Tikhonov & Arsenin 1977), realized by computing 3.5 where**

*Σ**α*is a small, positive, scalar (the regularization factor) and

**I**is the (

*N*

_{b}×

*N*

_{b}) identity matrix.

Implementation of equation (3.5) will, with sufficient regularization, always result in an answer (an estimated ‘bubble spectral generation rate’, ** Ψ**). The question is to what extent that answer is accurate, which cannot be ascertained without independent measurements, and certainly is a question which should be asked whenever the environmental conditions differ significantly from the forward model on which the physics is based (Leighton

*et al.*1998

*b*, 2002; Boyd & Varley 2001; Manasseh

*et al.*2001). While such independent measurements can be conducted to validate the technique (Leighton

*et al.*1996, 1997), in the field they may not be available, and so tests should be performed to examine whether the estimated bubble spectral generation rate can be discounted. A minimum test (which tests the process, not the model) is to determine whether the measured acoustic spectrum is reconstituted, when the estimated bubble spectral generation rate is inserted into the forward problem. A general rule-of-thumb for such inverse problems is to apply only enough regularization to make the solution turn from unstable to stable (for example, using the smallest value of

*α*for which none of the radius bins has a negative bubble generation rate within it). One should also note that the inversion inherently assumes a zero bubble generation rate for all bubble sizes with a natural frequency outside of the bandwidth of the detected signal. Iterative solutions (for example, by replacing the zero bubble counts outside of that bandwidth with smooth extrapolations of the population trend) can be tested for stability; the tendency for the number of large bubbles to decrease with increasing bubble size in many bubble populations supports the logic of this approach (Leighton

*et al.*2010).

The use of more than one appropriately positioned hydrophone is important, since the location of the sound source must be localized to a sufficient degree to remove the inherent ambiguity between the magnitude and range of the leak: a loud, distant leak might be indistinguishable from a nearer, quieter leak, unless more than one hydrophone is deployed, or other forms of waves (such as in the sediment or pipeline, or at the water/sediment interface or seabed), with different wave speeds, are used.

At the heart of the inversion is the assumption that the ambient noise measured by the hydrophone array is generated by a superposition of emissions from single bubbles for which equations (2.3) and (2.4) (or modified versions thereof) are appropriate models. Other sources of emission are possible in given frequency ranges (including both non-bubble ambient noise (Medwin 2005) and bubble emissions which do not conform to the single-bubble model (Loewen & Melville 1994; Leighton *et al.* 2002)), although human judgement is good at recognizing the superimposed sounds of single bubbles. A more objective test can be carried out by using the array to examine the sound emanating from other locations at increasing ranges from the leak (Leighton *et al.* 2005), since the acoustic element will decay with the inverse square of range from the leak, while other components of noise may not. The fact that the acoustic emissions from a leak can be identified, and their attenuation from the source of the leak to different ranges be modelled, can be combined with the ability of the array to sample the sound field at various ranges from the putative leak to iterate the analysed signal such that the bubble-generated component is identified for processing. The effective range of a hydrophone depends on the ambient noise, oceanic conditions, and size of the leak (which can vary greatly), and these variables are explored in §4.

## 4. Results and discussion

Time-series data from three methane seeps were used to construct the radiated sound field, assuming a forward model based on equations (2.5) and (2.6). These data were constructed using bubble generation rates obtained optically at sea (Sauter *et al.* 2006; Leifer & Culling 2010). The rates were used to compute the overall number of bubbles expected to occur during the measurement. These bubbles were generated at moments randomly assigned throughout the time series (according to a uniform distribution). Similarly, the radius of each bubble was randomly assigned with a weighted probability to shape a pre-defined bubble generation spectrum, and their time domain signatures being computed using equations (2.5) and (2.6). The spectrum of the resulting time-series was then estimated using standard non-parametric spectral estimation techniques (Kay & Marple 1981). Figure 1 shows these spectra computed using data derived from the at-sea measurements of Leifer & Culling (2010) for a ‘minor’ seep in the Coal Oil Point seep field at a depth of 45 m in the Santa Barbara Channel, California. The spectra are estimated in two ways: using a simulated time-series (dashed line) and directly from equation (3.1) (solid line). The discrepancies between these two curves arises because of the stochastic character of time-series generation model and the finite duration of the simulated signals, compared to the infinite averaging implied by the theoretical model. The spectra indicate the ‘source level’ in the standard manner (Urick 1983). That is to say, they characterize the far-field acoustic properties of the source, by indicating the sound pressure level of a hypothetical equivalent monopole source at some reference distance (1 m is used throughout as this reference distance).

It is important to bear in mind that these predictions assume that the release of bubbles into the water column generates acoustic emissions as predicted by equation (2.6). Their purpose here is to provide a baseline for the inversion, and at-sea measured spectra may differ because either the input bubble size distribution omits key bubbles (some large bubbles will be too infrequent, and some small bubbles will be invisible, for the video system of Leifer & Culling 2010, although low-pressure gas release from the seabed may involve coalescence and so not produce so many small bubbles as, say, rupture of a high pressure gas pipe; Leighton *et al.*, 1991) ; or the assumptions underpinning equation (2.6) may be invalid for certain forms of gas release. For some situations (e.g. bubbling of gas from a ruptured subsea gas pipeline) equation (2.6) may be valid, but for other circumstances it will not. If the gas injection occurs elsewhere (deep in the seabed) and the release of bubbles into the water column only represents the transport of two phase liquids from depths into the sea, then such a release will not conform to the model here whereby gas injection into the water column occurs at the seabed. In a similar vein, some bubble release mechanisms (such as bubble generation through chemical decomposition or electrolysis, or the exsolution of previously dissolved gas from water transported to shallower depths or warmer temperatures) can be far quieter than generation of the same bubbles through injection. No data on the possible acoustic emissions are given in the papers used here to infer data on seeps, and so equation (2.6) will be used in the first instance to demonstrate the method, knowing that for some circumstances it will be appropriate. Confidence however can be gained from the fact that the seep studied using acoustics and optics by Leifer & Tang (2007) generated observable acoustic emissions.

The inversions are all conducted for hydrophones situated 5 m from the centre of each seep, and based on the spectra constructed from the artificial time-series that a hydrophone at their range would receive. These spectra include both the emissions from the vent (i.e. as described by source level spectra such as figure 1, but modified by propagation out to the hydrophone positioned 5 m from the centre of the seep) and ambient noise. For want of specific noise spectral levels (NSLs), to illustrate the principle, the ambient noise is assumed to be white with a spectral level of 45 dB re 1 μPa^{2} Hz^{−1} unless otherwise stated (Urick 1983). Other studies have used separate hydrophones to estimate the ambient noise level (Leighton *et al*. 2005) and the ambient noise spectrum can be estimated in a robust fashion (Leung & White 1998). One can estimate the range, *r*_{d}, at which such a seep is just detectable. This is achieved by defining *I*_{s} to be the peak signal intensity in the frequency range of interest and assuming that detection occurs if this peak exceeds the noise spectrum in that frequency band (*I*_{n}) at least by an amount (Δ). Assuming that a spherical spreading loss of occurs between the range of 1 m (the reference position for the source spectral level) and range *r*_{d}, then the value of *r*_{d} for which detection occurs is given by
4.1
(in reality, topography and thermally or bubble-induced sound speed variations may produce spreading models other than spherical). Imposing a detection threshold of 6 dB for Δ results in a maximum detection range of 7 m for the minor seep considered here. Replacing the single omnidirectional hydrophone by an array of *N* hydrophones can increase this range. If it is assumed that the array is configured such that the ambient noise on the sensors in uncorrelated, but that the acoustic signal from the leak is perfectly correlated, then using simple delay and sum beamforming the detection range will be increased by a factor of √*N* (Urick 1983). One option is to employ a four hydrophone array, the elements of which could be arranged to remove left/right ambiguities. This would potentially double the detection range.

Figure 2 illustrates the result of applying the inversion to the estimated spectrum, based on 60 s of data gathered at a range of 5 m and embedded in a realistic level of ambient noise. The inversion is conducted using 100 radius bins in the range 0.5 mm to 10 mm. It is clear that the fit is good over a limited range of bubble sizes, but that errors occur for bubbles smaller and larger than these sizes. The ability to compare the result of the inversion with the actual bubble population used as input to create the spectra of figure 1 makes this apparent, but such information is not available in the field when monitoring real seeps. Crucially, local measurements can readily indicate the bubble size range over which the fit will be reliable, as follows. The frequency range where data from the hydrophone (at range 5 m) is not dominated by noise (which has a stated level of 45 dB re 1 μPa^{2} Hz^{−1}), is that frequency range where the source level exceeds dB. In figure 1 this occurs in the approximate range 1–5 kHz, which (from equation (2.5)) corresponds to a bubble size range of approximately 1.5–7.5 mm. Integration of the bubble generation rates across this ‘reliable’ range of bubble sizes can be used to estimate the gas flux from the inversion (a 0.75 compression factor is incorporated into the ideal gas law, to compensate for the compressibility of methane relative to an ideal gas, Sauter *et al.* 2006). Therefore, using only bubble data in the radius range 1.5–7.5 mm, the gas flux computed using inversion was 1.66×10^{−4} mol s^{−1}. This is accurate to within 4 per cent of the actual gas flux used as an input to form figure 1 (1.55×10^{−4} mol s^{−1}). The rarity of large bubbles, and the encapsulation of a relatively small proportion of the total gas volume in small bubbles for this population, contribute to this accuracy. This technique will now be applied to larger seeps, where the signal-to-noise ratio (SNR) at the hydrophone 5 m from the centre of the seep is greater. This might reasonably lead to the expectation of greater accuracy in assessing the bubble population, but this is not necessarily the case, since the dominance of bubble signal will probably generate a greater dynamic range in the hydrophone output, which makes regularization of the inversion matrix more difficult. In all these inversions, all frequency bands are included in the inversion and those associated with poor SNRs can be rejected after the inversion has been performed.

Figures 3 and 4 show the results of a similar exercise conducted on a ‘major’ seep measured by Leifer & Culling (2010) in the Coal Oil Point seep field, located at a depth of 45 m. The acoustic spectral levels generated by this larger seep are much greater (figure 3), peaking at 86 dB re 1 μPa^{2} Hz^{−1} at 1 m, giving from equation (4.1) a detection range of approximately 50√*N* m for an array of *N* hydrophones of the type described above. The bubble generation spectrum estimated from the inversion of 60 s of time series data is close to the original distribution (dashed line). There is still a tendency to over-count small bubbles because the SNR is lower at high frequencies for a hydrophone at 5 m range. In this case, the frequency range where the source level exceeds 59 dB re 1 μPa^{2} Hz^{−1} at 1 m, is approximately 1–9 kHz, corresponding to a bubble radius range of 0.8–7.5 mm. Note that in this case the input and estimated bubble size distributions (figure 4) stop at 5 mm, to be consistent with the original optical measurements of Leifer & Culling (2010). The allowance of a 7.5 mm upper limit is an artefact caused by energy emitted below the natural frequencies of the largest bubbles because of the finite bandwidth that is a consequence of their damping. In this instance, the estimated gas flux is computed from the inversion as 0.0023 mol s^{−1}, compared to the flux in the simulation which is equivalent to 0.0025 mol s^{−1} and so the inversion is accurate to within 11 per cent.

The final example considers the largest of the seeps for which we can infer bubble size data to use as input, located at the deep-sea submarine mud volcano (SMV) Håkon Mosby. This SMV releases gas over an area with a diameter of approximately 1500 m. In this work, it is the acoustic emissions from a single vent within the SMV which are considered. Sauter *et al.* (2006) do not give bubble generation spectra but from their photographic data it is possible to produce input for the scheme outlined in this paper. They consider a single small vent from a methane seep of 20 cm diameter at 1250 m depth, producing 1000 bubbles per second with a mean bubble radius of 2.6 mm. For want of measured bubble generation spectrum, a chi distribution with six degrees of freedom and scale parameter 2.45×10^{−4} is selected to realize a distribution with mean bubble radius of 2.6 mm. The chi distribution is selected because it ensures that the radii generated are always positive and provides a plausible uni-modal bubble size distribution. The spectrum of the synthesized acoustic signal, which includes ambient noise, is computed at an assumed range of 5 m and used as the basis of the inversion.

It is the loudest of the three seeps, peaking at 102 dB re 1 μPa^{2} Hz^{−1} at 1 m, such that from equation (4.1) its maximum detection range, using a single hydrophone, is 350 m. An even more substantial leak, which generates ten times the bubbles numbers modelled above with the same size distribution, would be detectable at a range of 3.5 km by a single hydrophone (although at such ranges the propagation model would need to include multipaths, topography, and absorption, and the thermal and bubble density gradients an SMV can generate). Replacing the single omnidirectional hydrophone by an array of *N* hydrophones, as described above, can increase these ranges to 350√*N* and 3500√*N* m. Since this is one seep of many in the volcano, in practice a hydrophone 5 m from the centre of this seep would probably suffer high levels of bubble-generated ‘noise’ from adjacent seeps in the volcano, and in the field the inversion would probably be done for all these seeps combined (i.e. the volcano as a whole) at ranges of km order. However, because bubble size data are only available for the one seep in the volcano, the above method will here be followed through for this one seep. The inversion is based on a spectrum estimated from 60 s of time series data and synthesized to replicate hydrophone data measured at a 5 m range. The resulting bubble generation spectrum is close to the input (figure 5, dashed line) for radii greater than approximately 1 mm. However, the inversion still over-estimates the generation rate for smaller bubbles. In this case, the acoustic energy exceeds the 59 dB re 1 μPa^{2} Hz^{−1} value for all frequencies above 800 Hz in the bandwidth considered (up to 50 kHz). The overall gas flux estimated via the acoustic inversion is 0.62 mol s^{−1}, compared to 0.69 mol s^{−1} which is the gas flux associated with the model, an error of 10 per cent.

In order to quantify the effect of noise level on the inversion results, simulations using the model of the major seep were conducted. The inversion is performed using different ambient NSLs. With the exception of the NSL the rest of the parameters of the model and inversion are the same as those used to compute figure 3. Figure 6 shows that the solution remains accurate for NSLs of up to approximately 60 dB re 1 μPa^{2} Hz^{−1}, corresponding to an ambient noise level at 5 kHz in sea state 6 in a deep water environment.

Finally, figure 7 shows the effects of incorporating some of the departures from the simple relations for bubble natural frequency (Minnaert 1933) and damping (Commander & Prosperetti 1989) which recent authors have raised as possibilities. For these tests, simulations of the minor leak were used with all other conditions matching those used to generate the results in figure 2. The time series data are generated from an input bubble population (shown in figures 2 and 6) assuming that nature has applied one of two departures from standard bubble theory: in the first case the bubble quality factors are halved from those used in standard bubble theory (after the suggestion of Walton *et al.* 2005), and in the second case the Minnaert relationship is adjusted (after Leifer & Tang 2007) such that a bubble of a given size has a natural frequency which exceeds by 20 per cent that predicted by Minnaert. The two time series produced are then inverted using the standard free-field theory for bubble natural frequency (Minnaert 1933) and damping (Commander & Prosperetti 1989), as if the experimenter seeking to invert the passive emissions from a seep were in ignorance that nature had applied slightly modified rules in generating that sound. Figure 7 demonstrates the estimated bubble generation spectrum and the true bubble spectrum and as expected additional errors are introduced as a consequence of the model mismatch. In calculating the gas flux, the same upper and lower bubble radius limits (1.5 and 7.5 mm respectively) are used as in figure 2, because the SNR is very similar.

When the quality factor is changed (dashed line in figure 7) after the suggestion of Walton *et al.* (2005), the bubble population is under-estimated, since the increased damping (reduced quality factor) results in a reduction in the acoustic energy emitted by each bubble, so that the overall power of the simulated time-series is reduced. As a result, in this instance, the gas flux is estimated (by integrating the results of the inversion over the radius range 1.5–7.5 mm) to be 9.52×10^{−5} mol s^{−1}, which is in error by almost 38 per cent from the input value of 1.55×10^{−4} mol s^{−1}.

The dotted line in figure 7 shows the result if the experimenter inverts using Minnaert's relationship (roughly [*f*/kHz]×[*R*_{0}/mm]≈5.3 at this depth), whereas nature generated the sound using roughly [*f*/kHz]×[*R*_{0}/mm]≈1.2×5.3 (after the suggestion of Leifer & Tang 2007). As expected the bubble radii estimated from the inversion are around 20 per cent smaller than the input (the peaks at around 2 and 3 mm in the input being reduced to around 1.6 and 2.4 mm, the peak at approx. 4.3 mm on the input having been smeared out). For these data, the gas flux estimated from inversion is 2.44×10^{−4} mol s^{−1} (an error of 60% relative to the model value of 1.55×10^{−4} mol s^{−1}).

An additional potential source of error can arise in the case of dense plumes (Deane 1997). If the attenuation through the cloud/plume is large, then the acoustic signature of a bubble located on the side of the plume distant to the sensor will be attenuated as it passes through the cloud, whereas bubbles located at the edge of the plume towards the sensor will not suffer attenuation through this mechanism. Assessment of the ability of this effect to cause undercounting of bubbles should be conducted, e.g. by calculating the attenuation caused by the population estimated by the inversion (or the one used as input if, as here, it is available). Here, the acoustic attenuation of a plane wave propagating through the plume with the greatest gas flux, i.e. the SMV seep, was computed (Commander & Prosperetti 1989). For this seep the acoustic attenuation through the plume was computed to be less than 0.25 dB. On this basis, the shielding effect of the plume has not been included in the model, or used to iterate the first estimation of the bubble population to generate a refined estimation. If this were required, the methods described by Deane (1997) and Deane & Stokes (2010) provide a method for incorporating this effect in situations where shielding is significant. Those authors found a more significant effect than that seen here because they considered the bubble population generated from above by breaking waves, where turbulence keeps many bubbles in the water column long after they cease to ring, leading to highly attenuating clouds which persist in the propagation path for long periods (Farmer & Vagle 1989; Medwin & Breitz 1989; Phelps & Leighton 1998). The low level of attenuation calculated for the seep plumes in this paper is in part owing to the fact that bubbles are generated from below and rise rapidly away from the base of the seeps where the bulk of the emissions are generated. At these frequencies, attenuation owing to particulate matter raised by the seep will similarly be negligible (Brown *et al.* 1998; Richards *et al.* 2003*a*,*b*).

## 5. Practical applications

The methodology described here to quantify gas flux requires that the range to the source of gas (leak/seep/plume) be estimated (equation (2.8)). Particular care needs to be taken to overcome the ambiguity between dense small clouds and disperse larger clouds (White *et al.* 2004), which is an acoustic equivalent of Olber's paradox. The mechanism by which such an estimate can be obtained is highly application dependent. There are very different practical constraints which can apply depending on the problem under consideration. To illustrate the possible solutions we detail three example applications: continuously monitoring a carbon sequestration site, continuous monitoring of a gas pipeline and studying a methane seepage site (figure 8).

The monitoring of carbon sequestration sites and gas pipelines requires a system which is deployable for extended periods, and capable of detecting when bubble generation occurs. Since plumes from both are hopefully infrequent indicators of failure in these man-made structures, passive acoustics can (if false alarm rates are low) provide a power-conserving trigger for alarms (with batteries potentially supported by energy harvesting), which then lead to more power-hungry confirmation systems (e.g. active sonar, optics, etc.). Passive acoustics arrays can provide first estimates of the plume location (and hence range to an individual hydrophone for quantification of gas flux) by estimating bearing using standard correlation methods (figure 8*a*,*b*). Once hydrophones had been placed in such a way as to gather passive acoustic data, if the low frequency bandwidth were sufficient the data could also be analysed for the low frequency collective oscillation, in the manner proposed by Maksimov (2004*a*,*b*, 2005), which would represent a complementary way of estimating gas flux using potentially the same hardware. Communication between units, and to a central control unit, is required, and the vulnerability, cost, and power and bandwidth constraints probably mean that coded information is sent in place of raw data.

The passive detection of waves with different intrinsic speeds (in pipe walls, water, seabed or piped gas) offers another route to locating the sound source. Simple active sonars can transmit above the pipe route and identify/locate leaks through perturbation in the amplitude, phase or speed of the transmission, or through a backscattered echo (figure 8*c*). Imaging active sonars can then confirm range and bearing from time-of-flight and beamforming (figure 8*d*) and, within limitation, quantify the gas flux (De Beukelaer *et al.* 2003; Sauter *et al.* 2006; Nikolovska & Schanze 2007; von Deimling *et al.* 2007; Westbrook *et al.* 2009), but the bubble size distribution would be an estimate limited by the bandwidth that is affordable (Maksimov & Sosedko 2002; Maksimov 2003*a*,*b*). However, the power required for an active capability can also be used to acquire bottom and sub-bottom sonar of the associated geological and gaseous features (Paull *et al.* 1995; Krastel *et al.* 2003; Naudts *et al.* 2008): such techniques, which can detect (and even quantify; Leighton & Robb 2008) gas release into the sediment before the gas reaches the water column, are particularly valuable for CCS facilities. Subsidiary information (e.g. pipeline route, probable failure points on structures) is also available. Some important leakages from CCS facilities may be less dramatic than gas pipeline ruptures, requiring different detection thresholds, and indeed other indicators (e.g. the effect of acidification on particular fauna and flora) may prove to be a useful monitor.

The issues are different for methane seeps, since these emit continuously and so an alarm based on passive acoustics is not usually useful. Methane seeps may physically be very large and many plumes are usually present simultaneously. Before fixed long-term monitoring installations are considered, it is probable that surveys will have quantified the leak, and sensor platforms will have ‘flown’ towards it to take data from various ranges (figure 8*e*,*f*). This would identify the range where SNR is suitable (but clipping avoided) for passive monitoring using the methods of this paper (figure 6), either to cross-validate other estimates of gas flux (from gas collection, active sonar, optical, etc.; figure 9), or to select the site for hydrophones for long-term change monitoring or averaging using passive acoustics (figure 8*d*), a low-power option to the active acoustic monitoring that has been achieved (Greinert 2008). In contrast to the mainly alarm-based role of passive acoustics for infrequent pipeline leaks, existing seeps must be characterized to provide input to model their environmental impact, which increasingly will require measures of the bubble spectral generation rates, bubble rise speeds (figure 9), hydrate and surface active species in the environment and on the bubble wall, etc. as inputs to models of plumes (Bettelini & Fanneløp 1993) and their acoustics (Maksimov & Sosedko 2009).

## 6. Conclusions

This paper describes a method by which probably the simplest, commonest and most robust of oceanic acoustical instruments, the hydrophone, can be used to provide quantitative measures of gas leakage from carbon storage facilities, methane seeps, and gas pipelines of a realistic size that they generate enough bubbles to be detected remotely. The equipment is suitable for long-term remote monitoring. The inversion has a single unknown parameter, *R*_{ε0i}/*R*_{0} (two unknowns in the unlikely event that simulations of the particular experiment show that the finite ring-up time is important) which can be evaluated using this theory from field trails, and refined with the developing theoretical base. However, the accuracy of the estimated gas flux depends on how well the circumstances of the leak match the assumed conditions. These are that the signal to be analysed (e.g. the hydrophone record after appropriate manipulation to subtract noise and—if necessary—identification of bubble signatures, Leighton *et al.* 1998*a*, etc.) conforms at source to the models of equation (2.3) or a modified version thereof, and that it is mediated by the environment (e.g. through modelled absorption, attenuation, scattering, diffraction and interaction with bathymetry) in a way that can be modelled. There are standard methods for such modelling (Jensen *et al.* 1994).

The inversions undertaken show an ability to estimate the bubble spectral generation rate and gas volume generation rates, although errors occur at the ends of range of the data, here most notably at the small bubble end of the spectrum. The estimated gas volume generation rates were within 11 per cent of the rates used as input for the assumed conditions (noise being particularly important). The louder the seep, the greater the detection range, here up to 700 m for a 4-hydrophone array, i.e. twice the 350 m detection range computed from the SMV seep based on a single sensor. The closer the observation is to the maximum detection range, the worse will be the SNR. This will produce a tendency for inaccuracy in the inversion, but other factors (such as difficulties in regularizing the matrix) might produce inaccuracies if the SNR is too great. Note that the largest seep considered here (the single SMV vent) is small compared to the usual flux from commercial pipeline leaks that are considered to be measureable, and indeed its detection and quantification (to within 10% accuracy in figure 5) is at least two orders of magnitude more sensitive than current model-based techniques for large, long pipelines (T. E. Bustnes 2011, personal communication; W. Postvoll 2011, personal communication).

The technique is only as valuable as the extent to which the physics in the model matches the physics at sea. Vulnerable assumptions, which will have to be tested at sea, cover both the bubble generation (in terms of to what extent the generation produces the assumed acoustic emission) and the propagation (in terms of absorption, topography, multipaths, etc.). Some studies have suggested departures in terms of measured damping and frequency, and this paper showed that if such effects existed at sea and went unrecognized, they could significantly decrease the accuracy of the inversion. However, insufficient work has been done to determine whether they will occur at sea, and further studies are required to separate confirmed effects and their sources from artefacts of the measurement. At sea, there could be a wealth of phenomena (such as the effect of the chemical constituents on the bubble damping and stiffness) not covered by theory, while in the laboratory there are often many phenomena (such as tank reverberation) which can affect these relationships (Leighton *et al.* 2002). Effects which are significant on single bubble dynamics (such as the generation of nonlinearities, Leighton *et al.* 2004) or damping (Ainslie & Leighton 2009) can, when incorporated into populations of bubbles to determine the overall effect, have only modest influences. The effect on gas flux estimates when the suggested perturbations to the Minnaert relationship and bubble damping were included in the forward problem (to generate the acoustic spectrum) but neglected in the inversion (as if the experimenters were in ignorance that such an effect had occurred) were assessed here and found to be large, indicating the importance of identifying real factors affecting the passive acoustic emissions from the ocean bubbles in question, and distinguishing them from artefacts (e.g. in the laboratory). Other forms of acoustic inversion (such as determining bubble populations from acoustic absorption, Leighton *et al.* 2004) might be more prone to changes in damping.

## Acknowledgements

The authors are grateful to Trond Erlend Bustnes of Statoil for advice and comment, and to Willy Postvoll of Gassco AS for information on the leak detection capability of current techniques used by industry.

- Received April 6, 2011.
- Accepted September 16, 2011.

- This journal is © 2011 The Royal Society