Broadband Acoustic Inversion for Gas Flux Quanti ﬁ cation — Application to a Methane Plume at Scanner Pockmark, Central North Sea

The release of greenhouse gases from both natural and man ‐ made sites has been identi ﬁ ed as a major cause of global climate change. Extensive work has addressed quantifying gas seeps in the terrestrial setting while little has been done to re ﬁ ne accurate methods for determining gas ﬂ ux emerging through the seabed into the water column. This paper investigates large ‐ scale methane seepage from the Scanner Pockmark in the North Sea with a new methodology that integrates data from both multibeam and single ‐ beam acoustics, with single ‐ beam data covering a bandwidth (3.5 to 200 kHz) far wider than that used in previous studies, to quantify the rate of gas release from the seabed into the water column. The multibeam data imaged a distinct fork ‐ shaped methane plume in the water column, the upper arm of which was consistently visible in the single ‐ beam data, while the lower arm was only intermittently visible. Using a novel acoustic inversion method, we determine the depth ‐ dependent gas bubble size distribution and the gas ﬂ ux for each plume arm. Our results show that the upper plume arm comprises bubbles with radii ranging from 1 to 15 mm, while the lower arm consists of smaller bubbles with radii ranging from 0.01 to 0.15 mm. We extrapolate from these estimates to calculate the gas ﬂ ux from the Scanner Pockmark as between 1.6 and 2.7 × 10 6 kg/year (272 to 456 L/min). This range was calculated by considering uncertainties together with Monte Carlo simulation. Our improved methodology allows more accurate quanti ﬁ cation of natural and anthropogenic gas plumes in the water column.

Active acoustics, specifically the use of multibeam echo sounders, has been commonly used for seep detection in the last decades (Greinert, Lewis, et al., 2010;Xu et al., 2020) and has been used to map both natural and anthropogenic ebullition sites worldwide (Greinert, Lewis, et al., 2010;Greinert et al., 2006;Leblond et al., 2014;Nikolovska et al., 2008;Ostrovsky, 2003;Ostrovsky et al., 2008;Urban et al., 2017;von Deimling et al., 2015;Westbrook et al., 2009).Echo sounders also have the advantage of being able to work in any body of water regardless of visibility unlike optical techniques.Gas seeps in sonar data commonly appear as readily identifiable medium/strong reflectors-within the water column, sometimes referred to as "gas flares".Using multibeam echo sounders, the position and shape of these flares can be mapped (Greinert, Lewis, et al., 2010;Urban et al., 2017).By mapping the shape of these flares, observing the angle they make with the seabed, and knowing the tidal velocity, one can predict the vertical velocity of the bubble cloud.There is a simple relationship between ascent velocity and bubble size, and hence, the dominant bubble size can be estimated (Toramaru, 1989).
In order to gain an estimate of gas flux via active acoustics, single-beam (single frequency) echo sounder data has been used (Bayrakci et al., 2014;Greinert, McGinnis, et al., 2010;Römer et al., 2014;Shakhova et al., 2015;Veloso et al., 2015;von Deimling et al., 2010).This is done by first modeling the theoretical return pulse strength from bubbles of different sizes based on the frequency of the acoustic source and the depth of water in the area.Then by observing the mean signal strength from within the plume, an estimate of bubble size distribution can be made.Crucially, this can only be done if the ambiguity is ignored, since when a given scattering strength is attributed to a bubble, there is always more than one bubble size that can scatter that frequency strongly (Leighton et al., 2004).Consequently there is an inherent ambiguity in the gas flux estimated by a technique that only uses data containing a single frequency.This ambiguity exists even when only a single bubble is being measured in free field (Leighton et al., 1996) and becomes much greater if there are many bubbles (as here) or the bubbles are contained within a structure (Baik et al., 2014;Leighton et al., 2012).From this distribution the flux of the plume can be estimated (carrying forward any inherent ambiguity) using the calculated rise speeds of bubbles (Greinert & Nützel, 2004;Greinert et al., 2006;Greinert, Lewis, et al., 2010;Leblond et al., 2014;Nikolovska & Schanze, 2007).Greinert et al. (2006) used both single-beam and multibeam data to estimate the dominant bubble size at different depths in the water column.This method has also been used to make observations of the temporal variations of plumes and their interaction with the thermocline (von Deimling et al., 2015).However, the modeling used in this method requires very accurate measurement of water column physical properties as well as bubble rise velocity.
In an attempt to establish a technique to directly quantify gas flux from active acoustic data, Greinert and Nützel (2004) demonstrated that (within the confines of a specific seep, constrained to remove the inherent ambiguity in the acoustic inversion) there is a direct relationship between the volume backscattering strength of a single-beam pulse and the flux rate of a seep, using a controlled release site and a horizontal acoustic array.However, this relationship varies with the dominant bubble size meaning it is site specific and must be re-established at every new seep via empirical measurements (Greinert & Nützel, 2004;Leblond et al., 2014).This approach was used by Nikolovska et al. (2008) in the Black Sea, using a Remotely Operated Underwater Vehicle (ROV) to collect physical flux measurements alongside a horizontally mounted sonar system, and by Bayrakci et al. (2014) in the Marmara Sea, using a rotating bubble detector (BOB) to reveal temporal variations in the gas flux of surrounding seeps.While this technique is appropriate for long-term measurements of single seep sites, it is intrinsically flawed for widespread quantification of multiple seeps as an empirical measurement of flux is required to make such an estimate.Furthermore, it assumes that conditions do not change (e.g., significantly larger bubbles are introduced through a new fracture in the sediment or infrastructure casing) in such a way as to make the gas flux quantification erroneous through the above-mentioned ambiguity.
The existence of the inherent ambiguity is therefore probably the most significant shortcoming of the existing acoustic techniques, that is, ebullition sites contain a bubble population with a wide range of radii (Veloso et al., 2015).Crucially, passive acoustic techniques do not contain inherent ambiguities in the acoustic inversion: They only contain acoustic uncertainties, which are less troublesome.To be specific, each bubble emits energy in a known frequency band relating to its size, depth, etc.The uncertainty in the amount of energy emitted by a given bubble being only a result of the paucity of data, which will be reduced as more data is taken (Leighton & White, 2011).In contrast, quantification of the gas flux by active sonar contains an inherent ambiguity, in that a given bubble can scatter strongly at resonance, and when it is also much larger than resonance (Leighton et al., 2004).As such, a single-frequency echo sounder can never unambiguously quantify the gas flux without additional measurements, for example, passive acoustics, optical methods, or gas collection using bottles, to remove the ambiguity.Measuring scattering across a range of frequencies, which does not cover the resonant frequencies of the bubbles present, leads to an ill-conditioned inversion problem; that is, the errors in the measurements are vastly magnified, leading to solutions that are unreliable.Physically relevant regularization of the solution is needed in order to provide usable solutions (Leighton et al., 1996).Furthermore, the above-mentioned active methods tend to rely upon scattering models for bubbles, which assume the bubble is small relative to the insonifying wavelength.For the size of bubble that we are looking at and the frequencies of most imaging sonars, this condition is not true.This leads to errors in two ways: First, the calculation of the damping associated with each bubble can be erroneous (Ainslie & Leighton, 2011, 2009); second, the assumed increase of scattering cross section with increasing bubble size (a trend that is only valid for bubbles larger than resonance only do as long as the bubble radius remains much smaller than an acoustic wavelength) breaks down (Salomatin & Yusupov, 2005;Thuraisingham, 1997).Accurate determination of the bubble population, and hence gas flux, can only be determined if the backscatter response is determined for all significant bubble sizes, and this requires the use of a broad range of acoustic frequencies.Typical radii of bubbles emitted from the seabed tend to be in the range of 1 to 15 mm  (Veloso et al., 2015) whose resonant frequencies are from 800 Hz to 12 kHz.While there is merit in using single-frequency imaging (at, e.g., 18 kHz;Xu et al., 2014) to identify the location of seep sites, single-frequency systems cannot determine the bubble population or the gas flux accurately.Even a multifrequency system that did not cover the range of bubble resonances (from below the resonant frequency of the largest bubble present, to above the frequency of the smallest bubble present) will contain inherent ambiguities, and if all the frequencies in a multibeam system are higher than the resonance of the larger bubble present (the convenient option given the frequencies in off-the-shelf multibeam sonars), then the equations in the simultaneous set mentioned above are not independent and cannot be solved to determine the variables (the number of the bubbles in each size bin) unambiguously (Berges et al., 2015;Leighton & White, 2011).Currently, researchers have been using optical methods for quantification of small plumes such as a single bubble stream (Veloso et al., 2015), but this is impractical for analyzing larger emission sites.Little work has been completed on quantifying the emissions from large methane plumes from active pockmarks, which may extend over a diameter of 200 m in the water column, or understanding the gas bubble upwelling process.
This paper proposes a solution to this dilemma using two elements.First, we employ a wider range of frequencies than previously used, so that more of the bubble resonances are encompassed.Second, it assumes a form for the bubble size distribution, further constraining the solution and effectively regularizing the inversion.We combine data from three sonar systems, spanning a wide frequency range, 2.5 to 200 kHz, to calculate methane flux at an actively venting pockmark in the North Sea.Scanner Pockmark complex comprises two large pockmarks (∼200-300 m in diameter, Figure 1), which are 15-20 m deep depressions in a relatively flat seabed in water depth of 150 m.Pockmarks are submarine gas escape structures commonly found in basins globally and often associated with active hydrocarbon systems.Despite first being observed in the 1960s, the variability and controls on gas emissions are poorly understood.The evolution of the resulting gas plumes in the water column is closely linked to the overall mechanism of gas leakage from pockmarks, making a greater understanding of plumes essential for better understanding natural seep sites.In order to determine the bubble size distribution of the gas plume and quantify the gas flux within it, we first use multibeam imaging to detect the plume structure and dimensions, then we present a volume scattering strength matching model utilizing iterations of bubble mean radii and standard deviation to match observed strength of single-beam data in the function of frequency ranging from 3.5 to 200 kHz for each depth.Next, a sea current modulation function is applied to integrate the instantaneous bubble rise velocity, estimated at the time of observation.Finally, we apply a depth-dependent number of bubbles and size distribution for methane gas to convert these volume flow rates to mass flow rates.

Data
The data in this survey were collected from the RRS James Cook during September 2017.Three hull-mounted sonar systems were employed: a Kongsberg EM710 multibeam echo sounder, a Kongsberg SBP120 sub-bottom profiler, and a Simrad EK60 single-beam echo sounder.The transceivers were orientated vertically downwards for the entire study.The EM710 multibeam echo sounder worked on frequency range 70 to 100 kHz with beam width of 1°.The SBP120 worked on a single-beam of wideband frequency (2.5-6.5 kHz) centered at ∼3.5 kHz with a beam width of 3°.The EK60 echo sounder transmits a single-beam of five different monochromatic frequencies: 18, 38, 70, 120, and 200 kHz, with beam widths of 11°, 7°, 7°, 7°, and 7°, respectively.The pulse length of the SBP120 was set to 40 ms, and the pulse length of the EK60 at frequencies 18, 38, 70, 120, and200 kHz were set to 2048, 1024, 512, 256, and256 μs, respectively.The multibeam data is used to observe the structure and dimensions of the plume, while the sub-bottom profiler and single-beam data are used to measure the acoustic scattering properties of the plume.The sampling rate of the SBP120 and the EK60 were set to 20.48 and 25 kHz, respectively, which makes it possible for the target strength calculation in a reverberating volume V (m 3 ) with 1 m vertical resolution.

EM710 Multibeam Data
The EK710 multibeam system imaged the methane plume from the western Scanner Pockmark.Figure 1shows the transects across the plume.By filtering out background water column noise, it was possible to extract the gas flare and recreate it as a 3-D model, recording its positional data, height, lateral extent, and width; example results of this process are shown in Figure 2 good agreement with tidal direction, the axis of which runs roughly north to south, as predicted by Cazenave's FVCOM model (Cazenave et al., 2016).The plume height varied between 39-145 m above the seabed while the lateral spread varied between 5 and 210 m.
Over the recent decades, numerous methane plumes in different ocean regions have been investigated and the occurrence of multiple arms has been noted on several occasions (Dissanayake et al., 2018;Gentz et al., 2014;Leifer et al., 2017;McGinnis et al., 2006;Ruppel, 2011;Sommer et al., 2015;von Deimling et al., 2015).Examination of Figure 2 reveals that the plume selected here exhibits a clear forked structure with two distinct arms.This is presumed to be a result of two dominant bubble sizes escaping from the pockmark, with the larger bubbles rise faster, creating the upper arm while the smaller bubbles rise more slowly creating the lower arm.
The multibeam data also allow us to map the surrounding seafloor topography, revealing the Scanner Pockmark as being 10-20 m deeper than the surrounding seabed, which is at a depth of roughly 150 m (Figure 1).It was not possible to clearly map the plume within the crater due to the increased reverberation, likely caused by internal reflections and active gas venting.

SBP120 and EK60 Single-Beam Data
Calibrated, single-beam data from the SBP120 and EK60 were collected along the transects A-D illustrated in Figure 1.38,70,120,and 200 kHz was collected along the four profiles A-D across the Scanner Pockmark plume using an EK60 system with built-in calibration.This was augmented by data collected from a 3.5 kHz (2.5-6.5 kHz) chirp sub-bottom profiler.An example of the plume imaged on one single-beam system is shown in Figure 3. Plume data were extracted by filtering out background water column noise data, based on the simultaneously collected multibeam data, leaving only the acoustic signal associated with the gas venting from the pockmark.Individual acoustic anomalies were removed if they were connected to the seafloor, or single, isolated, and vertically elongated stack of high acoustic energy above noise level.Additionally, the  1).The plume is orientated in the same direction as the tidal flow (i.e., in a north-south direction).The distinct forked shape of the plume can be observed.Plume lateral extent is colored from white at the base to black at the upper surface.and −30 dB; (a) 3.5, (b) 18, (c) 38, (d) 70, (e) 120, and (f ) 200 kHz.10.1029/2020JC016360

Data Processing, Modeling, and Flow Rate Estimation
This section describes how the observed target strength data is used to determine the depth-dependent bubble size distribution and gas flux of a plume (Figure 5).Acoustic detection and identification of gas plumes can be used to quantify the bubble flow rate if a number of acquisition parameters and assumptions about the physics of methane gas seepage at the seafloor and the surrounding environments are made (Veloso et al., 2015).The multibeam data are used to determine structure of plume arms and the corresponding dimension in depth.The single-beam target strength data are used to derive the observed volume scattering strength in depth.To quantify the bubble size distribution and gas flux, we develop an inversion algorithm, which iteratively matches the modeled and measured volume scattering data.For each depth of interest, the shape of the bubble size distribution is parametrized by a log-normal probability density function (PDF), with a further parameter defining the total number of bubbles.As mentioned in section 2, the Scanner Pockmark produces two dominant bubble sizes, and we incorporate this into our model.

Beam Data Processing
The intermediate frequencies of each data set are smoothed to create an observation of volume scattering strength as a function of depth and frequency.We denote the received target strength at frequency f of backscattering ping n as Ts n ( f ) (dB), then the volume scattering strength Vss r ( f ) (dB/m 3 ) can be expressed as (Johanneson & Mitson, 1983) where N p is the total number of scatterers in a fragment of volume, and the reverberating volume V is computed as where h i is the vertical height and S i is the scanning area of interest in the horizontal plane.Considering the propagation loss (PL( f )) in the acoustic channel, the volume scattering strength of gas bubbles Vss( f ) can be expressed as (Smailes, 1976) Essentially, the PL( f ) is the sum of two terms: the geometrical loss (PL g ( f )) and the absorption loss (PL α ( f )) (Li, White, Bull, & Leighton, 2019): Here we assume a spherical spreading model for the geometrical losses, and the absorption loss is calculated from Thorp's formula (Harris & Zorzi, 2007;Li, White, Bull, & Leighton, 2019;Ochi et al., 2008;Urick, 2013).
Taking into account the propagation loss, the volume scattering strength of gas bubbles as a function of frequency can be extracted.

Modeling
The model of the acoustic scattering from the bubble plume combines three basic components: (1) the model of the backscattering cross section of a single bubble, (2) an assumed shape of the bubble size distribution, and (3) a method to compute the volume scattering strength.Each of these three elements is detailed in a subsequent subsection.

The Acoustic Backscattering Cross Section of a Single Bubble
The backscattering cross section of a single bubble is relatively well established when wavelength of the ensonifying sound field is significantly greater than the bubble radius, that is, kr≪1 (sometimes referred to as the "long-wavelength" condition) where k is the wave number equal to 2πf/c w with f representing the frequency of the acoustic wave.When using many commercial imaging sonars to examine bubbles from seeps, this condition is frequently violated.For example, for the highest frequency used in this study, 200 kHz, then k is approximately 800 m −1 .To keep kr less than 5%, say, in order to make the "long-wavelength" formulations valid, the seeps should emit no bubbles larger than radius 60 μm for a 200 kHz beam.This maximum allowable bubble size, to keep the "long-wavelength" formulation valid, decreases with increasing frequency, for example, kr ¼ 0.05 for bubbles of 20 μm radius when f ¼ 600 kHz, by no means the highest frequency used to quantify gas from seeps.Given most measurements of seeps show bubble radii that are at least two orders of magnitude larger than this limit then the "long-wavelength limit" is not justifiable.
The gas flux from a seep is dominated by the gas carried in the largest bubbles, so to estimate such fluxes it is most important to accurately model the scattering from these large bubbles.As discussed in Ainslie andLeighton (2009, 2011), when the condition kr≪1 cannot be relied upon then one needs to take considerable care.In particular, it is necessary to ensure that expressions for the damping terms, arising through three mechanisms: acoustic radiation, viscous, and thermal damping, also do not rely upon assuming kr≪1.Further, the expressions for the cross sections need to be corrected from the prediction of the formulation of the long-wavelength limit (which erroneously predicts that the scattering cross section increases quadratically with increasing radius).They in fact approximately plateau (onto which resonances are superimposed), which is the prediction from detailed modeling (Thuraisingham, 1997).
The expression we shall consider for the backscattering cross section, σ bs , is (5) This is adapted from Ainslie andLeighton (2009, 2011) to include the final factor, which was proposed by Thuraisingham (1997).This expression implicitly includes radiation damping, with the effect of the other two damping mechanisms (viscous and thermal damping) being combined into a single damping factor, β 0 .This formulation provides a consistent approach to incorporating radiation damping into the backscattering model, something which, as Ainslie and Leighton (2011) showed, cannot be achieved using dimensionless damping coefficient, which is the prevailing approach (Veloso et al., 2015).In Equation 5, the frequency ω 0 is defined through the solution of the equation where R{•} denotes the real part of a complex number.Under specific circumstances (when the process is 10.1029/2020JC016360 Journal of Geophysical Research: Oceans isothermal or adiabatic), this frequency corresponds to the resonance frequency of the bubble; however, in general, this is not the case.The complex parameter Ω, seen in Equation 6, is defined through where ρ liq is the density of the liquid surrounding the bubble (kg/m 3 ), τ is the surface tension (N/m), and P gas is the pressure of the gas inside the bubble (Pa), which can be expressed as follows: where P atm is the atmospheric pressure (Pa), g is the acceleration due to gravity (m/s 2 ), p v is the vapor pressure for water, and d is the depth of the bubble.In Equation 7, Γ represents the complex polytropic index (Ainslie & Leighton, 2011) with γ representing the specific heat ratio, and the parameter X being defined as and the thermal diffusivity, D p , of the gas in the bubble can be expressed as in which the K gas is the thermal conductivity of the gas within the bubble and C p is the specific heat capacity of the gas at constant pressure.The density of the gas ρ gas can be computed using (Leighton, 1994) where R is the gas constant and T is the temperature (K) and M m is the molar mass of the gas.
The two remaining damped effects (thermal and viscous) are included in the model ( 7) through the combined damping coefficient β 0 defined as where β th and β vis are the thermal and viscous damping coefficient (s −1 ).Further, expressions for these two quantities can be obtained as follows (Ainslie & Leighton, 2009, 2011): where I{•} denotes the imaginary part of a complex number and η S is the shear viscosity of the liquid (Pa•s).The form for the viscous damping has been a matter of some discussion, with some Journal of Geophysical Research: Oceans authors favoring the inclusion of the effects of bulk viscosity (Love, 1978;Veloso et al., 2015); however, the later analysis of Baik (2013) highlighted flaws in the previous work and recommended the use of Equation 14b.
While this model captures much of the physics of acoustic scattering from bubbles in the large wavelength limit, it should not be regarded as complete.It still relies on the assumption that the bubbles are spherical, which for large bubbles will not hold true and can affect the backscattering cross section (Ostrovsky et al., 2008;Salomatin & Yusupov, 2005).Parameters used in the bubble backscattering cross section computation are summarized in Table 1.

Bubble Size Distribution Assumption
To estimate the bubble size distribution for each plume arm, a log-normal distribution (Johnson et al., 1994) is used as an appropriate bubble size distribution to match the plume bubbles (Veloso et al., 2015): and Thus, for each point at which the inversion is applied, we have three parameters to match: the mean radius r b in Equation 16, the standard deviation ς b in Equation 17, and the number of bubbles per unit volume N b .
The mean radius r b is related to the frequency f peak corresponding to the peak value of the volume scattering strength Vss peak ( f ) for each depth; the deviation ς b is related to the curvature C of volume target strength curve as a function of frequency f; and the number of bubbles N b is related to the amplitude of the volume scattering strength Vss( f ).The three parameters are initialized at the beginning of the iteration process.

Modeled Volume Scattering Strength
Assuming the backscattering of all bubbles at depth d are uncorrected, the modeled volume scattering strength d Vssð f Þ (dB) is the sum of the backscattering strength of the individual bubbles in radius bins centered on [r 1 , … , r end ], given by where N b (r n ) is the number of bubbles with radius r n per unit volume, following the bubble radius PDF p b (r) in Equation 15.For a series of frequencies f ¼ { f 1 , … , f end }, we obtain a vector of d Vssð f Þ.

Linear Inversion
One existing approach to quantifying gas from backscattered acoustic signals is based on linear inversion techniques.Such methods have been considered in cases when no bubbles are assumed to be resonant (Nikolovska et al., 2008) and without that restriction (Berges, 2015;Muyakshin & Sauter, 2010;Veloso et al., 2015).These methods are based on Equation 19, which can be expressed in matrix form: where the elements of the matrix A are the backscattering strengths, the column vector x contains the number of bubbles per unit volume within a size bin, and b contains the linear volume scattering strengths.Assuming the matrix A is of full rank and the number of radius bins is equal to the number of ensonifying frequencies, then this is a square system of equations with a unique solution, which can, in principle, be solved through matrix inversion.If the number of radius bins is less than the number of frequencies, a least square solution can be obtained, while if the number of radius bins exceeds the number of frequencies the problem is ill posed and an infinite number of solutions exist. 10.1029/2020JC016360

Journal of Geophysical Research: Oceans
One problem is that since only a small number of frequencies are typically used to ensonify a cloud, one can only estimate bubbles in a small number of radius bins.Further, the matrix A can become ill conditioned as the off-diagonal terms can become large compared to the diagonal elements.This is because both resonant bubbles and large bubbles generate high levels of scattering, so while the diagonal elements in A may be large, so too are the regions corresponding to large bubbles ensonified by high frequencies.This ill conditioning in A means that during the inversion process small errors are greatly magnified.This can be mitigated by imposing prior constraints on the problem, in the form of regularization and by ensuring that the solution is always nonnegative.
In this work we eschew the use of linear inversion and at the outset impose constraints on the assumed bubble size distribution, which leads to a nonlinear optimization problem for which cannot be solved within a linear framework.

Matching Procedure
Rather than adopting a least squares approach to minimize the difference between the observed and modeled volume scattering strengths we shall use a curve matching strategy.Such an approach allows one to match a curve across the frequency interval at a large number of points, rather than solving the problem at isolated points.There are several curve matching techniques that have been proposed, including the Smith-Waterman algorithm for sequence alignment (Gribskov & Robinson, 1996), the B-spline fusion technique (Xia & Liu, 2004), the Discrete Curve Evolution (Bai et al., 2007), and the optimal alignment method (Sebastian et al., 2003).Here we adopt a method based at the optimal alignment curve matching.
The iteration procedure for each plume arm is shown in Equation A1.We first identify the plume arm structure, measure the dimension for each identified arm in depth, compute the observed volume scattering strength Vss( f ) in depth, and prepare coefficients and environmental parameters collected in the experiment as shown in Table 1.For the matching process, we must initialize the bubble radius r b ð0Þ, the standard deviation ς b (0), and the total number of bubbles N b (0).The bubble radius r b ð0Þ is initialized from the plume upwelling velocity v v as described in Equation 21; the standard deviation ς b (0) is initialized as 1 mm; and the total number of bubbles per unit volume is initialized as a positive integer (here we use 100 for the upper arm possessing big bubbles, and 10,000 for the lower arm possessing small bubbles).The initial radii, r b ð0Þ, is selected to be 0.05 and 5 mm, for the lower and upper arms respectively.To accelerate the matching, one may need to adapt these initial values according to the observed target strength as a function of depth and frequency.
For each iteration n of the curve matching method, we calculate the volume scattering strength Vss( f ) as a function of frequency f through Equation 19.From the calculated Vss( f ) curve, we find frequency of the peak of the curve, f peak ðnÞ, the maximum absolute curvature, j Ĉmax ðnÞj, and the value of the peak volume scattering strength at that point, d Vss peak ð f ; nÞ.The magnitude of the difference between the modeled and observed volume scattering strength can be computed, that is, Δ d Vssð f ; nÞ ¼ j d Vssð f ; nÞ−Vssð f Þj.If the size of this difference is minimized (e.g., on average less than a threshold Th 1 ¼ 1 dB and the largest difference is less than a threshold Th 2 ¼ 3 dB) in a number of iteration steps (e.g., 50), then the iteration is stopped; otherwise, the parameters (r b , ς b , and N b ) are updated according to the recursions shown in Appendix A1.
After the iterative matching process, we obtain estimates of the mean radius r b , standard deviation ς b , and the number of bubbles N b as a function of depth in each of the plume arms.These three parameters define the PDF of the bubble size distribution as a function of depth, so that at any depth one can compute the gas volume and the gas flux.

Measurement 3.3.1. Measuring Plume Upwelling Velocity
In order to compute the gas flux, one needs to know the amount of gas at a given depth and the velocity of the gas.Individual bubbles rise through a liquid as a result of buoyancy, at a rate called the bubble rise velocity.A plume of bubbles also create motion of the surrounding water, creating a circulation (upwelling), this is the plume upwelling velocity and represents the velocity of bubbles in the plume, which is required in the flux calculation.
To estimate the plume upwelling velocity, we use the plume slope angle and modeled sea current speed.The average slope l p (highlighted in Figure 2) is obtained by measuring the height and extent of the plume.The 10.1029/2020JC016360

Journal of Geophysical Research: Oceans
slope of the plume varies with depth, tide, and current (Sündermann & Pohlmann, 2011).However, our multibeam data (Figure 2) suggests that the plumes observed here rise at an approximately constant angle and we use that angle to estimate a constant plume upwelling velocity.
We assume that the horizontal displacement of plume is entirely controlled by the current.Thus, we assume the relationship/slope angle θ between the horizontal displacement X h and vertical displacement X v of the plume is equal relationship between the horizontal velocity v h (the current) and the plume upwelling velocity v v .The plume slope is then given by Using Equation 21, the average plume upwelling velocity near the pockmark values ranging from 10 to 15 cm/s.These values correspond to the bubble rise velocities for bubbles with radii in the range 1-6 mm (Park et al., 2017).This is consistent with our choice of an initial mean bubble radius in the upper arm.

Gas Volume Estimation
The plume is assumed to have an ellipsoidal cross section in the horizontal plane as observed from the multibeam data (Figure 2).The major and minor axes of the ellipse are denoted D l and D s , which can be measured from the 3-D multibeam data.We consider the gas in the plume in terms of horizontal slices of constant height (here we use 1 m).The scattered signal measured at the single-beam echo sounder consists of contributions from a volume, which is approximately cylindrically shaped oriented along the axis of the beam shown.The length of the cylinder being c w τ/2 (where τ is the pulse duration) and the diameter of the cylinder is 2htanðβ=2Þ, where h is the depth and β is the beam width (Figure 6).Assuming a horizontal cross section of the plume is homogeneous, having the same properties as the observed beam fragment, we can apply our findings appropriately to represent the entire horizontal cross section of the plume.
Based on the estimated bubble size distribution, the gas volume b V p (L) within 1 m thick section through the plume can be approximated using Figure 6.Geometry for converting gas volume from reverberating volume area to plume volume in depth.The single-beam scanned area is a fragment of the plume cross section.The plume horizontal cross section is considered to ellipse as observed from the multibeam data shown in Figure 2. Note that the alongship width and athwartship width are function of the longest diameter and shortest diameter of the ellipse relying on the ship direction. 10.1029/2020JC016360

Journal of Geophysical Research: Oceans
Figure 6 shows the geometry for converting gas volume from reverberating volume area to plume volume in depth.The calculated gas volumes are only a fragment of the gas volume in the whole gas plume arm at their corresponding depth (or horizontal cross section).We measure the size of the horizontal area in the reverberating volume fragment V according to the beam width.With the measured plume dimension S h , and the gas volumes for each fragment b V p , we calculate the gas volume b V h for each horizontal cross section at each depth with h v ¼ 1 m thickness: where S p is the horizontal dimension of the volume fragment.

Gas Flux Determination
Because of the interaction between the plume arms and the sea current, the rise velocity of the small bubbles in the lower plume arm is forced to be similar to that of upper plume arm containing larger bubbles.The mean rise velocity of the bubbles in the upper subplume are calculated using Equation 21.The estimated plume gas volumes V g and the bubble rise velocity, we obtain the gas fluxes b f g (L/min) of the upper plume arm and the lower plume arm in depth: where h v is the volume thickness, which here is equal to 1 m, considering the pulse length and resolution.

Results
Primary observations of this data set (Figure 7) show strong volume scattering strength at 3.5-18 and 70-120 kHz.For larger bubbles that rise faster than smaller bubbles (due to increased buoyancy), the movement direction of bubbles is closer to vertical.The target strength of the plume is integrated as volume scattering strength in depth for 1 m thickness and smoothed as shown in Figure 7. Two bubble clouds are visible, one in 3.5-18 kHz and the other in 38-200 kHz, and with a somewhat blurred border between them.At the low frequencies (f< 18 kHz), the clouds are connected to each other without big gaps, while at high frequencies (f> 38 kHz), the clouds are more separated from each other compared to those at low frequencies.This is consistent with arm structures observed from the multibeam observations, with small bubbles (radii < 0.2 mm) producing the peak at around 120 kHz and large bubbles (radii > 0.2 mm) producing the peak at much lower frequencies (Figure 7, left column).
Using the model matching approach, section 3.2.5, we obtained the scattering profiles shown in Figure 7middle column.This process also yields estimates of the parameters defining the bubble size distribution as a function of depth.The difference between the modeled and observed volume scattering strength is shown in Figure 7 right column.For all these cases, we successfully matched the scattering strength with only small difference remaining.This process also yields estimates of the parameters defining the bubble size distribution as a function of depth.To verify the gas flux change in depth, we compare the results to the predictions from a numerical model, specifically the Methane Individual Bubble Impact (MIBI) model (Dewar, 2016).

Plume Structure Identification
The two-arm structure that we observe for the plume is consistent with that presented in the literature (Dissanayake et al., 2018;Gentz et al., 2014;Greinert et al., 2006;Leifer et al., 2017;McGinnis et al., 2006;Ruppel, 2011;Sommer et al., 2015;von Deimling et al., 2015).It is proposed that the observed plume structure is a consequence of two dominant peaks in the bubble size distribution.The plume was observed multiple times from different directions, and the two-arm structure is consistently observed (see Figure 2).Acoustic data available for volume scattering strength analysis are at water depths 39-73 and 86-145 m (Figure 7).

Bubble Size Distribution
Identifying the structure of plume is one of the important elements in quantifying gas flow rate.Another important issue is related to bubble size distribution of each plume arm.Here, we determine the bubble 10.1029/2020JC016360 Journal of Geophysical Research: Oceans Figure 8a shows the PDF of the upper and lower plume arms at depth 65 and 145 m, respectively.From the estimation, bubbles in the upper arm possess radii mainly between 1 and 15 mm, while bubbles in the lower arm possess radii mainly between 0.01 and 0.15 mm.In the upper arm, there are more bubbles at 145 m than that at 65 m at all radii.In the lower arm, there are more relatively large bubbles (0.025-0.15 mm) at 145 m than that at 65 m, while there are fewer relatively small bubbles (<0.025 mm) at 145 m than that at 65 m.

Gas Flow Rate Quantification
The Scanner Pockmark plume was imaged on multiple profiles in different directions and the volume extent of methane bubbles in the water column is well constrained by multibeam data.Using the water column volume mapped from the multibeam data, we could extrapolate the volume scattering strength-derived measurements from single-beam data from the four profiles.We accept that our transect-derived estimates of bubble density and distribution may be a simplification of the 3-D plume.
With the measured sizes and the maximum gas volumes for each fragment, we calculate the gas volumes b V h for each horizontal cross section at 1 m intervals in depth.With the estimated plume gas volumes, we obtain the dominant gas fluxes of both the upper and lower plume arms (Figure 9).From the calculation of gas flux for each size interval of bubbles, we obtain the relative gas flux contribution for each bubble size interval.It shows that the highest contribution of gas flux is from bubbles at radii of about 8 mm, and the contribution of gas flux from the lower plume arm can be omitted as shown in Figure 8b.
The results described in Figure 9 allow the estimation of in situ instantaneous flow rates in the water column, and for the upper plume this is 1.56 × 10 6 kg/year (294 L/min) at 145 m water depth, while for the same depth, the lower arm flow rate is 2.6 × 10 4 kg/year (4.9 L/min).In this form of depth-based estimation, the upper arm contributes 98% to the gas emission, whereas 2% are from the bubbles in the lower plume.The gas flux determination results suggest that the upper arm with large bubbles dominates the gas flux of the seabed released methane from the Scanner Pockmark.In addition to the flow rate estimates close to the pockmark, we also estimate gas flow rate in the water column associated with four different tidal heights.While at different tidal heights, comparing that of Profiles A (0.4 m) and C (0.2 m), for example, the variation of gas flux  9a) and can be up to 1.0 L/ min (6.0 × 10 3 kg/year at depth 165 m in the lower arm, Figure 9b).
The gas flux gradually decreases from 1.56 × 10 6 kg/year (294 L/min) at depth 145 m to 6.9 × 10 4 kg/year (40 L/min) in the upper arm at depth 40 m as a consequence of bubble dissolution.In the lower arm, no obvious trend in gas flow rate is visible due to the intermittent emission of smaller bubbles.Overall, for the western Scanner Pockmark plume, the instantaneous flow rate is estimated to 1.59 × 10 6 kg/year (299 L/min) at depth 145 m, and 2.4 × 10 6 kg/year (400 L/min) at depth 165 m extrapolating the results in Figure 9 downward to the base of the pockmark.This instantaneous value may not be wholly representative of average flow rate, as in this study we have not considered tidal or seasonal variability.

Gas Flux Verification of Upwelling Methane Plume Using Modeling
To verify the gas flux evolution along the plume, we apply modeling of the methane bubbles known as the MIBI model (Dewar, 2016).This is a modified version of a CO 2 bubble model developed by Chen et al. (2009), recreated in a study to compare impacts of CO 2 and methane in the water column (Dewar et al., 2013).To simulate methane, the dominant gas properties on the bubble dynamics and dissolution have been used.
The results obtained from the MIBI model are shown as lines in Figure 9a.The model suggests that the dominant bubble radii is somewhere around 4.5 mm (between 3.5 and 6.5 mm), as superimposed for comparison in Figure 9a.This size best matches the measurements made from Profile B. The dominant bubble radius in the plume appears to be around the 4.5-5 mm mark; this is slightly larger than the measured peak bubble radius of 3.5 mm.However, this is to be expected given that bubbles of radii up to 20 mm are measured in the plume.The MIBI model also predicts that the reduction in plume volume from dissolution and the bubble expansion from reduced pressure.The results at bubble radius between 4.5 and 6.5 mm match well with the acoustic measurements at Profiles A-D, validating the effectiveness of our approach.

Uncertainty Estimation and Discussion
To remove ambiguities, one must use frequencies both above and below those of the bubble resonances present.Most multibeam echo sounders have frequencies, which, unless one is looking at very deep seeps, are mostly higher than the bubble resonances present.To remove all ambiguities, the lowest frequency used must be lower than the resonance of the largest bubble present (calculated above to be around 1 kHz).With the  Journal of Geophysical Research: Oceans smallest bubbles presented calculated to have a resonance of 12 kHz, then with the energy available in this experiment from 2.5-6.5 kHz from the chirp, and 18 kHz from the echo sounder, we have obtained energy that is at frequencies that are less than the resonances of most of the bubbles present, but not below that of the largest bubbles present.This is an improvement, but does not reach the ideal of achieving a stiffness-controlled scattering from the largest bubbles present (Leighton, 1994).
When considering uncertainty in our calculation of total gas flux, we need to consider two components: the errors due to uncertainty in individual parameters, which are inputs into our calculation; and the propagation of this uncertainty in the model.We initially describe errors in individual parameters, before utilizing a Monte Carlo simulation approach to understand the uncertainty in the calculation of total gas flux.The Monte Carlo simulations are based on 1,000 repetitions, we define measures of uncertainty in the flow rate empirically for the pockmark by varying only one input parameter at a time, holding all others at constant values.
Table 2 describes the uncertainty in physical model parameters in the model.While the uncertainty in some physical parameters are small (e.g., temperature, sound speed in water, and salinity) and not significant, others are much larger.We focus here on discussion of those parameters that have significant uncertainty.The overall uncertainty of the flow rate estimation based on the applied physics is defined as a simple superposition (multiplication) of individual factors of uncertainty as follows.The temperature that we are using in the model is an averaged one of 8.1°C (or 281.25 K) in a range of measured temperature [7.6°C, 10.1°C] (or [280.75 K, 283.25 K]).The sound speed in the seawater was measured as between 1,483 and 1,489 m/s with an average of 1,485 m/s.For the shallow water scenario, we choose to calculate clean bubbles as we assume that gas hydrates are stable (Leifer et al., 2000).However, this clean bubble assumption in shallow water may not hold in all cases; thus, we include the dirty bubbles in the uncertainty estimation.Application of sea current and plume slope to determine plume upwelling velocity (then bubble rise velocity) remains a variation factor of 8% relative to the seabed.The matching difference between the modeled and observed volume scattering strength is limited within a threshold of 1 dB.
After the Monte Carlo simulation, we obtain the uncertainty of gas flow rate as follows.The temperature affects the viscosity and results in −2% to 0.4% uncertainty of the cross-section computation, with lower temperatures generally reducing the flow rate.The sound speed affects the wave number k and reradiation damping coefficient δ rad and can result in −0.4% to 0.3% uncertainty.Such uncertainty resulted from seafloor temperature, near seafloor salinity, sound speed, and seawater density in the shallow water shelf environments and impact on flow rate estimation was found to be nearly indiscernible.While the dirty bubble assumption reduces the flow rate and results in −21% lower gas flux estimation than that of the clean bubble assumption.The plume slope makes −11.5% to 12% uncertainty of plume upwelling velocity values, then the flow rate.A measurement of the overall uncertainty in the calculations can be defined by combining statistics of the range in estimated flow rate values and uncertainty from the theory of flow rate estimation.Totally, the cumulative uncertainty bounds on the average reported flow rates are −32% to 14%.We outline in the following our approach to define an overall uncertainty in the reported values of flow rates, summarized in Table 3.
Our estimated total instantaneous flow rates of 2.4 × 10 6 kg/year is a representative first-order value for the gas flow at the Scanner Pockmark in the central North Sea, and we propose a total uncertainty in the flow rate estimation of [−32%, 14%].However, if one assumes in the scattering model that kr≪1 (Thuraisingham, 1997;Veloso et al., 2015) then one estimates the flux as 1.3 × 10 6 kg/year and using the new model described here (section 3.2.1)that estimate becomes 2.4 × 10 6 kg/year.

Conclusions
In this paper, we developed a new methodology for calculating gas flux from a seabed seep using multibeam imaging, and quantification from single-beam echo sounders covering a broad bandwidth .We investigate a methane seep from the Scanner Pockmark in the North Sea and find that the plume in the Journal of Geophysical Research: Oceans water column is forked with two arms.The broadband methodology enables us to quantify gas flux with frequencies spanning the resonances of all the bubbles in the plume.It applies an iterative model to match the volume scattering strength from the water column for each of the plume arm.The matching results show that the upper arm comprises larger bubbles (1-15 mm in radius) and the lower arm comprises smaller bubbles (0.01-0.15 mm in radius).The total seabed methane gas flux is quantified to be between 1.6 and 2.7 × 10 6 kg/year (272 to 456 L/min) at the Scanner Pockmark.

Figure 1 .
Figure 1.Bathymetric map of the Scanner Pockmark complex with inset showing the position within the central North Sea.The position of four ship profiles (A-D) that acoustically image the methane gas plume within the western pockmark are shown.

Figure 2 .
Figure 2. The methane plume at the Western Scanner Pockmark imaged by the EM710 multibeam echo sounder (70-100 kHz) on four multibeam profiles (A-D; position shown in Figure1).The plume is orientated in the same direction as the tidal flow (i.e., in a north-south direction).The distinct forked shape of the plume can be observed.Plume lateral extent is colored from white at the base to black at the upper surface.
multibeam data allowed us to cross validate the position of plumes and ensure that the relevant target was being examined.

Figure 4
Figure4shows examples of the target strength collected by the single-beam systems.Each sonar data set consists of the target strength of the plume at a range of depth in response to frequencies from 3.5 to 200 kHz.

Figure 3 .
Figure 3. Direct 18 kHz single-beam observation from snap shots of plume for Profile C. With the single-beam data we are unable to observe the forking in the multibeam data shown in Figure 2 due to the 2-D profile orientation.

Figure 5 .
Figure 5. Flow chart describing the processing steps used to determine bubble size and gas flux from the input acoustic data.Blue blocks describe the processing of observed multibeam data, single-beam data, and tidal information; purple blocks describe the iterative volume scattering strength matching model; and green blocks describe the quantification stage.

Figure 7 .
Figure 7. Left column: observed volume scattering strength of gas bubbles as a function of depth and frequency for the four profiles across the Scanner Pockmark methane.Data input was the volume scattering strengths observed at frequencies of 3.5, 18, 38, 70, 120 and 200 kHz; intermediate values are smoothed from the available data.Middle column: matched volume scattering strength as a function of depth and frequency for the four profiles.Right column: difference between the matched and the observed volume scattering strength as a function of depth and frequency for the four profiles.After sufficient iterations, the mean and maximum differences between the matched and observed volume scattering strength for most of these profiles are limited to 1 and 3 dB, respectively.

Figure
Figure (a) Bubble radius distribution estimated from volume scattering strength matching at depths of 65 and 145 m, respectively.In the upper arm, bubble radii are predominantly in the interval (1 mm, 15 mm); in the lower arm, bubble radii are predominantly in the interval (0.01 mm, 0.15 mm).(b) Relative gas flux comparison for each bubble radius bin at depth 145 m.The highest gas flux is contributed by bubbles with radii of about 8 mm.

Figure 9 .
Figure 9.Estimated gas flux for each plume arm.(a) Upper plume arm; lines with markers on are measured gas flux for each profile; lines without markers on are from the MIBI modeling (Dewar, 2016) of the methane bubbles, which match the measured gas flux quite well in depth; 0.5-6.5 mm are bubble radius.(b) lower plume arm.

Table 1
Parameters Used in the Cross Section Computation

Table 2
Measured and Applied Parameters in the Estimation Approach

Table 3
Uncertainty Estimation of Gas Flux F g (L/min) Using Monte Carlo Approach

Journal of Geophysical Research: Oceans
Vss peak ( f ) peak value of observed volume scattering strength |C peak | maximum absolute value of curvature of observed volume scattering strength curve f peak frequency corresponding to peak value of observed volume scattering strength V ss modeled target strength in the function of frequency f V ss peak peak value of modeled volume scattering strength j Ĉpeak j maximum absolute value of curvature of modeled volume scattering strength curve f peak frequency corresponding to peak value of modeled volume scattering strength Δ V ss difference between observed and modeled volume scattering strength Δ V ss mean mean value of the volume scattering strength difference Δ V ss max maximum value of the volume scattering strength difference V p estimated gas volume of observed fragment V h estimated gas volume of an entire horizontal cross section with 1 m thickness F g estimated gas flux 10.1029/2020JC016360