Sonar gas flux estimation by bubble insonification : application to methane bubble flux from seep areas in the outer Laptev Sea

Sonar surveys provide an effective mechanism for mapping seabed methane flux emissions, with Arctic submerged permafrost seepage having great potential to significantly affect climate. We created in situ engineered bubble plumes from 40 m depth with fluxes spanning 0.019 to 1.1 L s−1 to derive the in situ calibration curve (Q(σ )). These nonlinear curves related flux (Q) to sonar return (σ ) for a multibeam echosounder (MBES) and a single-beam echosounder (SBES) for a range of depths. The analysis demonstrated significant multiple bubble acoustic scattering – precluding the use of a theoretical approach to deriveQ(σ) from the product of the bubble σ (r) and the bubble size distribution where r is bubble radius. The bubble plume σ occurrence probability distribution function (9(σ)) with respect to Q found 9(σ) for weak σ well described by a power law that likely correlated with small-bubble dispersion and was strongly depth dependent. 9(σ) for strong σ was largely depth independent, consistent with bubble plume behavior where large bubbles in a plume remain in a focused core. 9(σ) was bimodal for all but the weakest plumes. Q(σ) was applied to sonar observations of natural arctic Laptev Sea seepage after accounting for volumetric change with numerical bubble plume simulations. Simulations addressed different depths and gases between calibration and seep plumes. Total mass fluxes (Qm) were 5.56, 42.73, and 4.88 mmol s−1 for MBES data with good to reasonable agreement (4–37 %) between the SBES and MBES systems. The seepage flux occurrence probability distribution function (9(Q)) was bimodal, with weak 9(Q) in each seep area well described by a power law, suggesting primarily minor bubble plumes. The seepage-mapped spatial patterns suggested subsurface geologic control attributing methane fluxes to the current state of subsea permafrost.


Arctic methane
On a century timescale, methane (CH 4 ) is the next most important anthropogenic greenhouse gas after carbon dioxide (CO 2 ) (Forster et al., 2007).However, on a decadal timescale comparable to its atmospheric lifetime, CH 4 is more important to the atmospheric radiative balance than CO 2 (Forster et al., 2007;Fig. 2.21).After nearly stabilizing, atmospheric CH 4 concentrations are increasing again, although the underlying reasons remain poorly understood (Nisbet et al., 2015).Despite likely increasing future natural emissions from global warming feedbacks (Rigby et al., 2008) and anthropogenic activities (Kirschke et al., 2013;Wunch et al., 2009), many source estimates have large uncertainties with greater uncertainty in future trends.This is particularly relevant for Arctic sources where global warming is the strongest (Graversen et al., 2008).
Arctic continental shelf sediment accumulates five times faster than in other world oceans.Sedimentation for the Siberian Arctic Shelf where the six Great Siberian Rivers outflow has deposited organic carbon that approximately equals accumulations over the entire pelagic area of the world's oceans.This leads to the thickest (up to 20 km) and most extensive sedimentary basin in the world (Gramberg et al., 1983).
Sonar is the most common survey approach and has been used on concentrated seep areas covering ∼ 1000 m 2 in the North Sea (Schneider von Deimling et al., 2007, 2010;Wilson et al., 2015), far more dispersed and weaker seepage in the Black Sea of ∼ 2500 plumes in an area of ∼ 20 km 2 (Greinert et al., 2010), and offshore Svalbard where a few hundred plumes were observed in an area of ∼ 15 km 2 (Veloso et al., 2015).Sonar can also be used from remotely operated vehicles for the deep sea, e.g., Muyakshin and Sauter (2010) for the Haakon Mosby mud volcano (3 plumes).
Sonar has also mapped significantly larger and stronger seepage in the Coal Oil Point (COP) marine hydrocarbon seep field offshore California.The COP seep field covers ∼ 3 km 2 of active seabed in an 18 km 2 area releasing 10 5 m 3 CH 4 per day (Hornafius et al., 1999), likely composed of many tens of thousands of plumes.A single survey requires two days (Leifer et al., 2010).
ESAS seepage is on a dramatically larger scale with ∼ 30 000 plumes manually identified in just two transects (Shakhova et al., 2014;Stubbs, 2010).Seepage densities up to ∼ 3000 seep bubble plumes per km 2 were found transecting a single hotspot.Based on the hotspot size (18 400 km 2 ), an order of magnitude estimate suggests 60 million seep plumes for the hotspot alone.Two sonar survey transits of the ESAS required a month.

Study motivation
Given the ESAS seepage extent there is a critical need for new approaches to effectively, rapidly, and quantitatively survey seepage areas.Video is inadequate to survey extensive or widely dispersed seepage, a task for which sonar (active acoustics) excels.This study demonstrates an improved approach to quantify seabed seepage using in situ calibrated sonar-derived bubble fluxes.Bubble plumes were observed in the ESAS and offshore California.In combination, the in situ studies covered a broad range of flows and included finedepth resolution of near-source (growth) plume processes (California) and coarser resolution of plume processes to tens of meters where the plume is self-similar.
Both multibeam echosounder (MBES) and single-beam echosounder (SBES) data were collected in the ESAS, just MBES data of rising engineered bubble plumes were col-lected in California.These data were collected both as a depth-dependent calibration and to investigate the effect of multiple acoustic scattering in bubble plumes.
The calibration was applied to quantify in situ sonar observations of three natural seepage areas in the ESAS.Because the calibration bubble plumes and seep bubble plumes were different gases and from different source depths, bubble dissolution rates are different -i.e., for the same seabed mean volume flux, the depth-window-averaged volume fluxes are different.We demonstrate a first correction attempt based on a numerical bubble plume model between the calibration and natural seepage bubble flows.
The ESAS is the world's largest and shallowest shelf (covering 2.1 × 10 6 km 2 ) containing the largest area of submerged permafrost by far (Shakhova et al., 2010a, b).The ESAS is a seaward extension of the Siberian tundra that was flooded during the Holocene transgression 7-15 kyr ago (Romanovskii et al., 2005).The ESAS comprises ∼ 25 % of the Arctic continental shelf and contains over 80 % of global subsea permafrost and shallow hydrate deposits, estimated at ∼ 1400 Gt carbon (Shakhova et al., 2010a).This reservoir includes hydrate deposits of ∼ 540 Gt of CH 4 with an additional two-thirds (∼ 360 Gt) trapped below as free gas (Soloviev et al., 1987).

Permafrost degradation
ESAS subsea permafrost degradation allows the release of sequestered CH 4 to the shallow ocean and atmosphere and has been changing in response to glacial and interglacial Arctic warming (∼ 7 • C) and warming from the overlying seawater (∼ 10 • C) since inundation in the early Holocene with additional ESAS seawater warming in recent decades (Biastoch et al., 2011;Semiletov et al., 2013Semiletov et al., , 2012;;Shakhova et al., 2014).
Recent observations of offshore permafrost (Shakhova et al., 2014) show that the ESAS permafrost lid is perforated, with year-round CH 4 emissions to the atmosphere coming from the sedimentary reservoir (Shakhova et al., 2010a(Shakhova et al., , 2015)).This migration feeds a vast marine seep field entirely in shallow waters where emissions contribute directly to the atmospheric budget (Shakhova et al., 2014).
There are important geologic controls on the subsea permafrost's thermal state.On millennia timescales, increasing temperatures of the overlying bottom seawater affect subsea permafrost through heat transfer and salinization (Shakhova et al., 2014(Shakhova et al., , 2015;;Soloviev et al., 1987).Geologic control also arises from heat transport by large Siberian rivers, which drives bottom water warming and is proposed to control the distribution of open taliks in coastal ESAS waters (Shakhova et al., 2014).Global warming enhances terrestrial riverine heat including from ecosystem responses, degradation of terrestrial permafrost, and increased river runoff.Warm riverine runoff drives a downward heat flux to shelf sediments and subsea permafrost (Shakhova and Semiletov, 2007;Shakhova et al., 2014).
Subsea permafrost degradation is greatest in the outer shelf waters (deeper than 50 m) where submergence occurred first, such as the outer Laptev Sea where models predict mostly degraded permafrost (Bauch et al., 2001).Riverine input to the Laptev Sea also supports the formation and growth of subsea thaw lakes and taliks, which are effective gas migration pathways to the seabed (Hölemann et al., 2011;Nicolsky and Shakhova, 2010;Shakhova and Semiletov, 2007;Shakhova et al., 2014).
Active seafloor spreading in the Laptev Sea, which is undergoing continental rifting, leads to strong geologic heat flow (85-117 m W m −2 ).The outer Laptev Sea is one of the few places where active oceanic spreading approaches a continental margin (Drachev et al., 2003) and correlates with the "hot" area crossed by the Ust' Lensky Rift and Khatanga-Lomonosov Fracture (Drachev et al., 2003;Nicolsky et al., 2012).Evidence of rifting is provided by hydrothermal fauna remnants documented around grabens (dropped blocks between faults) in the up-slope area that typically occurs along oceanic divergent axes (Drachev et al., 2003).Grabens in the ESAS are often linear structures that correlate spatially with paleo-river valleys.

Marine seepage fate and bubble processes
Marine seepage is a global phenomenon where CH 4 and other trace gases escape the seabed as bubbles that rise to-wards the sea surface (Judd and Hovland, 2007).These bubbles dissolve and deposit CH 4 in the water column, transporting their remaining contents to the sea surface -if they reach it (Leifer and Patro, 2002).The fate of these bubbles and their gas depends strongly on depth, size (Leifer and Patro, 2002), flow strength -plume synergies that include the upwelling flow velocity (V Z ; Leifer et al., 2009), and bubble surface properties like contamination (Leifer and Patro, 2002).
The fate of dissolved seep CH 4 depends most strongly on its deposition depth (Leifer and Patro, 2002) with CH 4 below the winter wave mixed layer (WWML) largely oxidized microbially (Rehder et al., 1999).In the shallow Coal Oil Point seep field, most of the CH 4 reaches the atmosphere directly (Clark et al., 2005) from mixing in the near field (Clark et al., 2000) and in the far (down-current) field when winds strengthen as typically occurs diurnally for coastal California.The same is true for the shallow ESAS where virtually all seabed CH 4 (dissolved and gaseous) is emitted in the WWML and escapes to the atmosphere directly by bubbles or through air-sea gas exchange by frequent storms (Shakhova et al., 2014).Even for deep-sea seepage (to ∼ 1 km), field studies show seep bubble plume CH 4 transport to the upper water column and atmosphere (MacDonald, 2011;Solomon et al., 2009) from plume processes and hydrate skin effects (Rehder et al., 2009;Warzinski et al., 2014).
Bubble size is important with most seep bubbles in a narrow range.Based on a review of 39 bubble plume size distributions (the most comprehensive published dataset to date), Leifer (2010) found that the vast majority of reported seep bubble plumes could be classified in two categories termed major and minor, with minor being the most common -see also the studies reviewed in Leifer (2010).Bubble plume size distributions ( (r e )) for minor bubble plumes are well described as Gaussian and largely composed of bubbles in a narrow size range, 1000 < r e < 4000 µm, where r e is the equivalent spherical radius.Major bubble plumes generally escape from higher-flow vents with a power law size distribution (Leifer and Culling, 2010).Most major bubble plumes are small, down to r e < 100 µm; however, plume volume is primarily found in the largest bubbles, up to r e ∼ 1 cm (Leifer et al., 2015b).

Sonar seep bubble observations
Seep (r e ) have been measured by video (Leifer, 2010;Römer et al., 2012;Sahling et al., 2009;Sauter et al., 2006) and passive acoustics.Passive acoustic (r e ) measurement has only been demonstrated for low-flow bubble plumes where the individual bubble acoustic signatures can be identified (Leifer and Tang, 2006;Maksimov et al., 2016).
Sonar interpretation is highly challenging, even for qualitative assessments of relative emission strength.For SBESs, there is geometric uncertainty (Leifer et al., 2010) et al., 2011;Wilson et al., 2015) where the sonar returns along multiple pathways, creating ghosts, shadow noise, off-beam returns, scattering loss, and other artifacts (Wilson et al., 2015).Note that if bubble spatial densities are sufficiently high for artifacts to occur between plumes, then the bubbles inside the plumes will produce artifacts inside the plumes.The vessel acoustic environment can be noisy and signal loss from scattering can also occur from suspended sediment and biota, often in layers.
There are many challenges to the quantitative derivation of bubble emission flux from sonar return, which at its basis relates to the interaction of sound with a bubble.For a single spherical bubble, the relationship has long been known with resonance given by the Minnaert (1933) equation: where f o is the resonance (or Minnaert) frequency, γ is the polytropic coefficient, P is internal bubble gas pressure, and ρ is pressure.This is the frequency a bubble emits at formation.For nonspherical bubbles (r e > 150 µm) there is an eccentricity correction based on the bubble-axis-wave front angle.Bubble eccentricities vary from 1.0 for spherical bubbles to 2 or greater for r e > 3500 µm (Clift et al., 1978).For a single spherical bubble, the back-scattering cross section (σ ) near f o is (Weber et al., 2008) where f is the frequency and δ is the damping term, approximated as δ ∼ 0.03f 0.3 with f in kHz.For bubbles larger than resonance, σ varies within 5 or so dB; for bubbles smaller than resonance, it decreases precipitously -35 dB for a factor of 2 decrease in r e (Weber et al., 2014).Integrating over the bubble emission size distribution, which is the number of bubbles in a r e bin, passing through the measurement plane combined with the bubble vertical velocity (V Z (r e )), which is a function of r e over the measurement volume yields the total plume cross section if bubbles are acoustically noninterfering (no multiple scattering) and the bubble-sonar interaction is isotropic -i.e., σ B is independent of angle despite bubble eccentricity.However, bubbles are often clustered in close proximity in seep bubble plumes, which allows multiple scattering between bubbles that decreases σ significantly (Weber, 2008).Acoustic modeling of bubble clusters is complicated even for small spherical bubbles -e.g., see Weber (2008).Specifically, σ depends on (r e , x, y, z) in the plume, which is asymmetric with respect to currents and evolves as the bubble plume rises.Acoustic propagation across a plume varies with azimuthal angle because bubbles are compressible, leading to different bubble-bubble acoustic interactions.Bubble eccentricity also contributes an azimuthal angle dependency.Whereas artifacts, like ghosting, between plumes (not side lobe return) can be spatially segregated, this is not feasible for such artifacts inside the plume.Here the bubbles are within a few centimeters (∼ 10-20 r e ) of each other, such as bubbles in the wake of a larger bubble (Tsuchiya et al., 1996) and near the seabed, acoustic coupling leads to frequency shifts (Leifer and Tang, 2006) that decrease σ .

Field study areas
This study reports on the use of in situ engineered plumes for calibration of σ to derive quantitative flux rates using a MBES which was deployed in the Coal Oil Point seep field, offshore California in the northern Santa Barbara Channel, in the Kara Sea of the ESAS.We present the small fraction of collected Kara Sea and ESAS data that were cleared for publication.

Coal Oil Point seep field
A precursor study was conducted in the COP seep field (Fig. 1) prior to the Arctic field experiment to demonstrate time-resolved, 3-D seep monitoring by a scanning MBES (Fig. S1 in the Supplement).The rotator-lander was de-ployed ∼ 15 m from the center of Shane Seep, which covers an area of ∼ 10 4 m 2 in ∼ 20 m water depth and comprises on the order of 1000 individual vents or bubble plumes (Fig. 1b).The lander included a MBES (Delta T, Imagenex, Vancouver, Canada) and compass (Ocean Server, MA) on an underwater rotator (Sidus Solutions, CA) with an azimuthal rotation angle range of up to 270 • .The sonar produced a 260 kHz vertically-oriented 128-beam fan spanning 120 • tilted upwards to reduce seabed backscatter.Two in situ calibration air bubble flows were deployed ∼ 8 m from the lander at azimuthal angles beyond the active seepage area and were traversed during each sonar rotation cycle.Two rotameters measured regulated airflows from an onboard compressor to these two bubble plumes.

Arctic field campaign
Field data were obtained during an expedition onboard the R/V Viktor Buynitsky from 2 September to 3 October 2012 (Figs. 2 and 3).The Vicktor Buynitsky sailed from Murmansk to the Laptev Sea and the adjacent portion of the ESAS.The expedition's overarching goal was to improve the understanding of the current scale of ESAS CH 4 emissions in order to develop a conceptual model of CH 4 propagation from the seabed to the atmosphere, including assessing source strengths and their dynamics.
The calibration experiments were conducted in a region of no natural seepage and almost flat seafloor in the Kara Sea (Fig. 3) to reduce or eliminate off-beam acoustic seabed scattering.Water depth was 45 m and weather was favorable: calm sea with a wind speed of 1-3 m s −1 and wave height of 0.2-0.5 m with no significant waves (0 to 1 ball).Column profile temperature and salinity data were measured using a conductivity, temperature, and depth (CTD) profiler (SBE19+, Seabird, USA).Weather for the seep sonar survey was typical (3-4 storm events with wind speed > 10 m s −1 ).
The vessel was anchored during the engineered bubble plume experiments.Engineered bubble plumes were made from nitrogen supplied by a pressure tank on the vessel foredeck.A 70 m long, 12 mm diameter, 6 mm wall thickness gas supply tube was attached by a Kevlar rope to a heavy metal weight (∼ 30 kg) that ballasted against buoyancy of gas in the tubing and drag from currents.The supply tube was deployed to a 40 m depth (Fig. S3) and the rising bubble plume was observed with a MBES and SBES.The sonars were located near each other so that their beam coverage overlapped with the center beam focused on the end of the bubble stream.Bubbles were produced from a 4 mm diameter copper nozzle attached at the end of the gas supply tube.
Gas flow was controlled using standard flow meters.One port was connected to a PVC tube and a second port was connected to a two-way valve.The third port was connected to the gas tank through the gas manifold.The manifold consisted of a high-pressure sensor for the tank pressure and a low-pressure sensor for the outgoing pressure (5.5 bar).We used temperature-compensated differential-pressure sensors with a manufacturer-specified range of ±1 psi (equivalent to ±70 cm of water).The sensor has a manufacturer-specified accuracy and stability of ±0.5 % FSD (full scale deflection over the operating-pressure range of the sensor over one year, between 0 and 50 • C) and repeatability errors of ±0.25 % FSD.For the study, the gas flow was varied from 0.5 to 150 L min −1 at 5.5 bar (equal to the bubble outlet hydrostatic pressure).For each experiment, the gas flow was allowed to stabilize and then sonar data were recorded for ∼ 10 min.
The same MBES was used in the ESAS and COP seep field.The SBES was a SIMRAD EK15 SW 1.0.0 echosounder (http://www.simrad.com) at 200 kHz, with a 1 ms pulse duration at 10 Hz, a 26 • beam width, and a built-in calibration system.Sonar data, including seep bubble plumes, were recorded at an average survey speed of 4-6 knots.Sonar backscatter was calibrated using acoustic targets (SIMRAD, Denmark).Initial data visualization and processing used EchoView and Sonar 5 software (SIMRAD) for the EK15 echosounder.
Bubbles have high density contrast with water and thus are strong sonar targets that can be distinguished easily from the background (Fig. 4b).For the engineered bubble plume experiments, the wave-mixed layer (WML) extended to ∼ 35 m depth with upper water ∼ 3.5 • C warmer than deeper water (Fig. 4a).
Sonar data analysis and visualization was performed with custom MATLAB routines (MathWorks, MA) that first georectified each ping and then assembled the data for each experimental run into a three-dimensional array of depth (z), transverse distance (x), and along track distance (y) or time (t) if stationary.

Seep and engineered bubble plume modeling
A numerical bubble propagation model was used to explore the relative dissolution rates for seep versus calibration bubble plumes and to calculate a volumetric correction factor that accounted for the difference.The bubble model is described elsewhere (Leifer et al., 2006;Leifer and Patro, 2002;Leifer et al., 2015b;Rehder et al., 2009).The model solves the coupled differential equations describing bubble molar content (Eq.3), size (Eq.4), pressure, and rise for each bubble size class in a bubble plume.These equations describe how sonar observations of bubble volume (size) relate to bubble mass (molar content).
Bubble dissolution or gas flux (F i ) for each gas species (i) is the change in bubble molar content (n i ) with time (t) driven by the concentration difference ( C i ) between the bubble and the surrounding water: where k B is the individual bubble gas transfer rate and depends on the gas diffusivity and r e , A is the bubble surface   area, H is the Henry's Law equilibrium, and P is the bubble partial pressure.Seep gases, such as CH 4 , largely outflow (positive F ) while air gases inflow (negative F ). F i depends on depth and r e through k B and also A (Leifer and Patro, 2002).Deeper bubbles of the same r e contain greater mass, allowing for longer survival and rise.Seep bubbles are seldom isolated (Leifer, 2010); thus, plume processes are important, including the upwelling flow which depends on the total plume volume flux (Leifer et al., 2009;Leifer et al., 2006).Another plume process is enhanced aqueous concentrations relative to the surrounding water, which enhances bubble survival (Leifer et al., 2006): were R is the universal gas constant, T is temperature, and n is the molar sum of all gases.This first term describes how the flux changes the bubble molar content and hence the change in bubble size with time (t).The second term describes how changes in hydrostatic pressure as the bubble rises (i.e., as depth (z) decreases) affect bubble size and depend on water density (ρ W ) and gravity (g).The denominator also includes the effect of surface tension (α) on pressure; higher pressure implies a smaller bubble.Unfortunately, bubble size distributions were not measured, and thus a typical minor bubble size distribution from the literature was used.Implications of these simplifying assumptions are discussed in Sect.4.4.
The model was initialized with a typical (Leifer, 2010) minor (Fig. 5a) for either methane or nitrogen bubbles, dissolved air gases at equilibrium in the water column, the observed CTD profile (Fig. 5b), and a 10 cm s −1 upwelling flow.V Z is an average value that is too low for the highest calibration flow and too high for the lowest (Leifer, 2010).

Engineered bubble plumes
Sonar return for the two (high and low) calibration plumes (Fig. S2) was thresholded above background (bubble-free water) and integrated for each beam during rotation across each calibration plume.The thresholded σ in a depth window then was fit with a linear polynomial of the log of the integrated sonar return over the plume versus height (h).As the bubble plume rose, σ increased -i.e., σ (h) was not constant (Fig. 6).Note, the change in volume for air bubbles over such short rise heights is negligible.This is evidence of decreasing bubble-bubble acoustic interaction as the bubbles rise and spread from turbulence (acoustic interactions decrease towards zero as the inter-bubble distances increase).Note that these data were uncalibrated and cannot be directly compared to the Arctic calibration data; this is presented to show the depth trends.
There is significant geometric uncertainty in SBES data, which is evident in the overlap in time of sonar returns for the calibration bubble plume (Fig. S4).This overlap results from current advection of the plume orthogonal to the page.MBES addresses this SBES deficiency.For example, the SBES loses the bubble plumes once they have risen into the WML, where currents often shift, but the MBES continues to observe them to 13 m depth, slightly below the vessel's draft.
The most common sonar return ping element is noise, which was isolated from the bubble plume signal by setting a threshold from the sonar return probability distribution function ( (σ )) at approximately −80 dB (Fig. 7a).The (σ ) weaker than −70 dB is clearly distinct from the stronger, but less common (lower ), bubble (σ ).Based on inspection of (σ ), a noise threshold value of −70 dB was selected (Fig. 7a, arrow), which provided a 5-8 dB transition between noise and bubbles.Obvious sonar artifacts, which can exhibit strong sonar return signatures, were masked by spatial segregation.Specifically, the plume center was identified at each depth and then filtered to ensure continuity with depth.Then, only samples within a specified horizontal distance from the plume centerline that tightly constrained the plume above the noise threshold were incorporated into the analysis.
For the engineered bubble plume experiments, plumes with volume flux (Q) from 0.019 to 1.1 L s −1 were created and observed by both SBES and MBES systems (Fig. 7).The contribution of bubble plume weak and strong sonar returns were investigated by their signature in (σ ).Specifically, (σ ) was modeled by a piecewise least-squares linear regression analysis of (σ ) = aσ (z) b .This model was then compared to expected trends in plume evolution of a rising bubble plume.Fit parameters are shown in Table S1 in the Supplement.Example data and fits for the 0.8 L s −1 plume shown in Fig. 9d-f for three depth windows (all below the WML).
For low-versus high-flow plumes, (σ ) was distinctly different, whereas (σ ) for the intermediate-flow plume exhibited characteristics of both low and high flows.A weak σ represents small bubbles, whereas a strong σ may reflect large bubbles or dense aggregations of small and/or large bubbles.As a bubble plume rises, the relative importance of small bubbles should increase as small bubbles disperse, spreading the weak sonar return over a larger volume.
(σ ) at the deepest depth for the weakest bubble plume exhibited a clear, two-part power law (Fig. 7c; Table S1) and remained constant as the bubble plume rose for the first 10 m, abruptly steepening in the next 5 m.This emphasizes the importance of smaller bubbles (b = −8, −7, and −12 for weak σ for the 45-40, 40-35, and 35-30 m depth windows, respectively).For the weaker bubble plumes (0.042 and 0.019 L s −1 , Fig. 7b and c, respectively), the strongest σ disappear completely at the shallowest depth, consistent with bubble plume dispersion, bubble dissolution, and strong currents.
(σ ) is bimodal for the deepest depth window for the highest-flow plume (Fig. 7d) with stronger σ more common relative to weaker σ than in the low-flow plume (Fig. 7c) and more common than "predicted" by extrapolating the weak σ power law fit (σ −10.7 ) to stronger σ (Fig. 7d and f, respectively).As this plume rose, (σ ) for weak σ decreased in relative importance while (σ ) for stronger σ remained constant.The power law exponent for the intermediate depth (b = −7.4) was less steep than for deeper (−10.7) and shallower (−8.4) depths.Thus, most of the evolution of (σ ) is from spatial expansion of weaker σ , i.e., smaller bubbles, while the denser, stronger σ bubbles remain relatively uniformly constrained with depth.The overall increase in σ with rise exhibits the same depth evolution as observed in the precursor COP study (Fig. 6) which featured strong plumes comparable to the strong plumes in Fig. 7d-f.(σ ) for the intermediate-flow plume (Fig. 7b) shares characteristics of both the high-and low-flow plume (σ ) -bimodal at the deepest depth with a pronounced strong σ peak in (σ ), like the high-flow plume, and evolving into a dual power law as the plume rises, as for the low-flow plume (σ ).Thus, (σ ) for the intermediate-flow plume evolved as it rose through the patterns of the strong to the pattern of the weak flow plumes.
These are point source plumes that disperse as they rise, thus bubble-bubble multiple scattering should decrease with height.With the exception of the strongest plume, plume rise decreases σ ; however, for the strongest flow plume, rise initially increases σ , similar to the behavior in the precursor study (Fig. 6) which was for comparably high flows albeit over depths much closer to the source.See Figs.S5 and S6 for sample MBES data for these flows.
The depth-dependent calibration curves (σ (Q, z)) were derived to account for the depth-evolving bubble-bubble acoustic interactions as the bubbles rose (Fig. 8).Specifically, σ above the noise threshold in spatially segregated boxes in each depth window is averaged over 7 min of sonar data for each flow to derive σ (Q, z).The MBES and SBES calibration datasets show saturation at high flow, similar to Greinert and Nützel (2004), which is evidence of bubble-bubble multiple scattering as shown in simulations by Weber (2008).For the high-flow cases, this likely includes sonar shadowing of more distant bubbles by nearer bubbles (decreasing σ ).At low flow, σ increases with increasing Q far faster than the linear addition of the number of bubbles.For example, for a flow doubling (Q = 0.02 to 0.04 L min −1 ), σ should increase by 20 log 10 (2), or 6 dB, yet increases are much larger.
σ (Q, z) shows a depth dependency in σ for both SBES and MBES systems (Fig. 8).For low-flow plumes, σ decreases with rise and is nonlinear with Q.In contrast, for high flows, both SBES and MBES systems saturate or are near saturation although there is significantly more variability in the MBES data.Saturation occurs when increased Q has a minimal to no increase in σ .Close inspection of the high-flow plume MBES data revealed undulations, which may have led to depth aliasing of σ in the 5 m depth windows.Although the high-flow calibration plumes are relevant for major seep bubble plumes such as in COP seep field (Leifer, 2010), the ESAS plumes studied were weaker.Thus, the strong calibration plumes are not discussed further.In contrast, the lowflow calibration plumes are comparable to typical minor bubble plumes (Leifer, 2010) and span the observed range of natural seepage in the study area.

Bubble dissolution rates and volume flux
These in situ σ (Q, z) were derived for application to seep bubble-sonar survey data and account for the bubble vertical velocity from buoyancy and upwelling flow.However, application of σ (Q, z) should account for the depth difference between the seep study area and the calibration plumes The Cryosphere, 11, 1333Cryosphere, 11, -1350Cryosphere, 11, , 2017 www.the-cryosphere.net/11/1333/2017/S2.
(70 versus 40 m) and different gas composition (seep gas was primarily CH 4 , while the calibration gas was nitrogen).Both these factors have non-negligible implications for the bubble dissolution rates of the two different plumes.These differences cause different bubble plume evolution and thus different volume height profiles.A volumetric correction factor was developed based on the ratio of the volume height profiles between a calibration and a seep bubble plume (same bubble size distribution) based on numerical bubble propagation model simulations (Fig. 9).
The numerical simulations show that for the first three 5 m depth windows, the depth-averaged total bubble plume volume (< Q z >) increases by 4.7, 15, and 29 % (Fig. 9b).This occurs primarily from decreasing hydrostatic pres-   sure and secondarily from oxygen inflow, while it shrinks from nitrogen outflow.Growth indicates that the balance favors hydrostatic over nitrogen outflow.
The size distribution of a minor seep bubble plume changes dramatically as it rises from a 70 m depth, with the smallest bubbles dissolving and the largest bubbles growing (Fig. 9d).Overall, air uptake and decreasing hydrostatic pressure largely balance dissolution for the first 50 m of bubble rise and < Q z > remains roughly stable (Fig. 9e); Q decreases by 0.7, 0.2, and 0.0 % in the first three 5 m depth windows.Note that stable Q does not imply constant total CH 4 bubble content, which continually outflows the rising bubble.
Combining the volumes from the two simulations provides the volume correction factors, 0.948, 0.868, and 0.775 for the 65-70, 60-65, and 55-60 m depth windows, respectively.Thus, the calibration plume Q averaged over the 35-40 m depth window is ∼ 5 % greater than the seep bubble plume Q for the 70-65 m depth window.

Natural seepage
The calibration function (σ (Q, z)) was applied to MBES and SBES Laptev Sea sonar data under strong-current conditions.Flux for the seep areas (Fig. 10) was mapped by averaging the seepage flux in the 65-70 m depth window in 1 m 2 quadrats after the application of the calibration function and the volume correction factor.The deepest depth window was chosen to better preserve the seabed location of emissions for spatial analysis.Three seep areas were surveyed, two weak and one strong, and all with numerous plumes.The MBES data illustrates the additional spatial information missing in SBES systems.For example, Seep area 1 in the SBES data (Fig. 10b) appears to show extensive diffuse seepage, which the MBES data (Fig. 10a) reveal arises from many low-flow discrete bubble plumes.
Seep area 2 was stronger than the other seep areas by an order of magnitude and clearly showed a northeast-southwest trend, which is apparent in all seep areas.Some of the striation patterns, primarily of the weaker returns, are consistent with the very strong currents detraining small bubbles out of the plume in the direction of the sonar beam fan.On a second, east-west leg, Seep area 1 was surveyed with currents not aligned with the sonar beam fan and does not exhibit striation.Further evidence of the effect of currents is shown in the sonar ping data (Fig. 10b vs. 10c and d), where Seep area 1 does not show the extreme tilt across beams as sonar data for Seep areas 2 and 3. Thus, the linear seep trends reflect geological control.
The mass flux (Q m ) occurrence probability distribution function ( (Q m )) was calculated for each seep area and showed that Seep area 2 contained the largest number of strong seep plumes followed by Seep area 3 and then Seep area 1 (Fig. 12).For the seep areas, (Q m ) for weak emissions asymptotically approached ∼ 0.1 mmol m −2 s −1 (2.5 cm 3 s −1 ) -the noise level.Thus, calibration flows (Fig. 8) were bracketed from the MBES data from the noise floor to the largest observed seep plume.
Seep area 2 exhibits both greater fluxes and a shallower power law than other seep areas (Fig. 12c).Furthermore, all seep areas exhibited positive anomalies or peaks in (Q m ) for stronger flux seepage.These peaks signify a preferred emission mode -i.e., multiple seeps with similar emission fluxes.For weaker seeps with good signal to noise ratios (Q m > 0.15 mmol m −2 s −1 ), the power law fits are nearly identical: 6.65, 6.27, and 6.80 for Seep areas 1, 2, and 3, respectively (Table 2).
Total flux in each seep area was determined by area integration and was 5.56, 42.73, and 4.88 mmol s −1 for the MBES data (Table 2).SBES-derived emissions were biased lower compared to MBES, by 3.7 to 36 % for the seep areas, with best agreement for Seep area 2.

Bubble-bubble acoustic interaction
We presented results of an in situ engineered bubble plume experiment to investigate the evolution of bubble plume sonar return for flows spanning two orders of magnitude.Table 1.Integrated depth-windowed methane flux estimates.This range was comparable from typical low-flow minor plumes to very strong high-flow major plumes (Leifer, 2010).Calibration plume sonar return increased strongly and nonlinearly with flux, ∼ 15 dB for a flow doubling from 0.02 to 0.04 L s −1 .This increase is faster than the 6 dB increase that would be expected by simply summing the sonar cross sections of the doubled number of bubbles.Instead, the increase suggests strong bubble-bubble acoustical interactions.Specifically, with increased flow, overall plume dimensions expand more quickly, leading to less bubble shadowing and shallower sonar occurrence probability distribution function slopes at the same height above the nozzle (Fig. 7).
In contrast to the overall plume dimensions (which includes smaller, more-dispersed bubbles), the dense core of large bubbles tends not to disperse and is largely insensitive to height (Fig. 7).Thus, for the dense core, increased flux increases bubble shadowing (multiple scattering) such that the signal of the additional bubbles is "blocked" (i.e., diminished) by other bubbles and sonar return becomes nearly independent of flow, i.e., saturated (Fig. 8a and b).Similar saturation is apparent in the data presented in Greinert and Nützel (2004) for an air bubble plume in far shallower water.Thus, the calibration data provides strong evidence of nonnegligible bubble-bubble acoustical interaction at both lowand high-flow rates.Furthermore, the relationship's nonlinearity is shown in the trend of σ (z, Q) as the bubble plume rises and disperses.Thus, bubble multiple scattering is significant even after the plume has risen 15 m.
Seep As a high-flow bubble plume rises, the weak σ portion of the plume representing small bubbles disperses, leading to an increase in the integrated σ as was observed in the Coal Oil Point and ESAS engineered plume data.In the COP seep field study, calibration flows extended from comparable to far higher flows than those in the ESAS and documented that σ (z) increased with height on fine-depth scales (Fig. 6).This was interpreted as due to decreasing bubble "shadowing" of more distant bubbles as the plume expands and becomes more diffuse during the plume growth or acceleration phase (Leifer et al., 2015a).As the ESAS calibration plumes rose, the sonar occurrence probability distribution function showed a strong influence from small-bubble dispersion as the plume expanded and an increase in the integrated σ (Fig. 8) As low-flow calibration plumes rise and disperse, σ decreases.Overlapping intermediate depth windows were evaluated and confirmed that this was not an artifact of plume oscillatory motions aliasing the return signal across the depth windows.The decrease in integrated σ with rise is (by definition) a decrease in scattered sonar energy, i.e., greater energy scatters back to the sonar when the plume is spatially denser.This could arise from a decrease in shadowing from scattering or dissolution; however, the bubble model showed that minor plume dissolution did not change overall plume volume significantly (Fig. 9), unlike the significant changes in integrated σ (e.g., Fig. 8c).

Bubble detrainment and bubble-bubble acoustic interaction
The artifact striations in the natural seep sonar data from currents are consistent with a non-negligible bubble-bubble acoustic interaction (Fig. 11).Specifically, seep bubble plumes were imaged for high currents that advected small bubbles out of the plumes into the down-current water.When detrained bubbles were in the beam fan orientation, they were observed, but not when the beam fan was perpendicular to the currents.For co-orientation, scattered acoustic energy interacts with nearby down-current bubbles, which remain in the beam.This arises because the cross-track beam dimension is very broad (120 • ) while the along-track beam width is very narrow (a few degrees).Thus, when cross-oriented, the sonar beam fan fails to image detrained bubbles.This provides clear evidence of bubble-bubble scattering at larger distances than the plume dimensions.

Bubble size distribution
Bubble size distributions have been reported for other ESAS seep sites (Shakhova et al., 2015), but the equipment to make bubble measurements was unavailable for this study.Bubble modeling was used to address the effect of the evolving bubble size distribution with flow in the application of calibration air or nitrogen (preferred for safety reasons over methane) bubble plumes to seep bubble plumes (Fig. 9).Thus, we applied a first approximation using a typical minor bubble plume size distribution.Clearly initializing the model with measured plumes would improve the accuracy of the volume correction factor and hence the sonar-derived flux.Still, the primary goal in our study is to demonstrate with a simple approximation that bubble size evolution matters and should not be neglected.
Although the simulations were conducted to correct between a nitrogen calibration plume and pure methane seep bubbles, if the seep bubbles contained other gases at nontrace levels, their outgassing could significantly impact bubble size evolution.In particular, CO 2 , which is far more soluble than CH 4 , can lead to rapid bubble size change, primarily in the deepest depth windows, e.g., see CO 2 plume simulation in Leifer et al. (2015b).Additionally, greater sensitivity arises from bubble plume depth (Leifer and Patro, 2002).Thus, the depth discrepancy between calibration and seep plumes should be minimized.Future calibration studies should also account for size distribution and upwelling flow with respect to flow rate.

Field comparison of MBES with SBES
The MBES and SBES systems were calibrated with the same nitrogen gas bubble plumes; thus, the two systems should agree in terms of flux observations.Calibration flows spanned very weak flows (Q = 0.19 L s −1 ) to very strong flows (Q = 1.1 L s −1 ).The low-flow calibration bubble plume (Fig. 8) was less than the seep field noise floor of the MBES system (Fig. 12a).In contrast, the high-flow calibration bubble plume was more than an order of magnitude greater than field observations.Field observations showed far better agreement between systems for Seep area 2 than the other seep areas (Table 2).This most likely relates to the greater relative importance of stronger seeps that are well above the noise level relative to the other seep areas.The calibration flows (Fig. 8) showed weaker sonar return for the SBES than for the MBES for the same flow.Geometric uncertainty likely played a role in the SBES negative bias.

Seepage spatial characterization
The seepage spatial map in the ESAS (Fig. 11) shares similarities with spatial patterns in the COP seep field (Fig. 1).Subsurface geologic structures control the seepage spatial flux distribution by creating the pathways through which seepage migrates to the seabed and ocean; seepage areas must occur where geologic structures allow.In the COP seep field, strong seepage areas are located at intersecting noncompressional faults and fractures (Leifer et al., 2010).Furthermore, these faults and/or fractures themselves are preferred migration pathways that connect subsurface reservoirs to the seabed, with seepage tending to manifest along their trend.
Two spatial trends were manifested in the ESAS seepage map (Fig. 11): one northeast-southwest of individual vents and the second a north-south elongation in Seep area 2. Both trends were aligned with the two weaker seepage areas.Furthermore, the northeast-southwest trend is apparent within Seep area 2. Here, fractures in submerged permafrost could play a similar role to the role of fault intersections in the COP seep field; however, more extensive seep area mapping is needed for validation and/or penetration sonar data that can image near-surface rock strata.On smaller length scales, there is an evident striation pattern in vent locations suggesting a subsurface linear geological control on meter-length scales.
High-flow seepage requires high permeability migration pathways, while low-flow seepage occurs along low permeability migration pathways if the driving pressure between the deeper reservoir and the seabed is constant across the active seepage area (Leifer and Boles, 2005).Thus, the stronger, more numerous, and extensive seepage emissions from Seep area 2 indicate higher subsurface permeability and subsurface connectivity with more numerous migration pathways than the other seep areas (Fig. 11).Seepage connectivity can be envisioned topologically as an inverted branched structure (Leifer and Boles, 2005) where central stronger seepage is surrounded (generally) by weaker seepage (Fig. S7).Given that permeability is inversely related to resistance in the migration pathways, stronger seepage is fed by migration along pathways with lower resistance (higher permeability), while weaker seepage is fed by migration along pathways with stronger resistance (lower permeability).The balance between seepage emissions for different migration pathways with a range of permeability underlies the flux probability distribution function (Fig. 12).
The seepage emissions map demonstrates similar geologic spatio-flux control.Specifically, weak seepage exhibited a b = −6.5 power law (Fig. 12), which describes the distribution between high and low permeability migration pathways.This argues that the shallow seabed structure (fracture, porosity, etc.) related to low permeability migration pathways is common across the areas, with the main controlling factor being the number of bubbles escaping per second per unit area of seabed.Note that although b is affected by bubble detrainment into the beam fan for Seep areas 2 and 3, Seep area 1 does not exhibit this effect yet has a similar b.
This power law does not extend to the largest seep fluxes, which manifest as peaks in the flux probability distribution function.Thus, higher-flow plumes could represent normal seabed structure failure (that governs the weak seepage) from stresses and/or talik melting, leading to focused high-flow migration pathways that help define where the seep areas lie.
In the Arctic, subsea permafrost degradation from heating both below (geologic, strongest in faulted zones) and above (riverine inputs and overall Arctic Ocean warming) creates migration pathways that manifest as seep spatio-flux distributions.The presence of active seepage in this region likely relates to these heat flows, with the hotspots likely related to taliks and/or subsea thaw lakes whose locations are controlled by linear geologic structures.In the ESAS, grabens are often linear structures, which often are correlated with paleo-river valleys, and could cause co-aligned fractures controlling seepage along linear trends.The similarity in the emission probability distribution power laws between seep areas indicates that subsurface permeability exhibits similar fractal distribution between the three areas.This argues for a similar formation mechanism, i.e., taliks.In this case, at the intersection of the two linear trends, fluid migration and thus heat flow are likely higher, leading to more rapid talik development providing high permeability migration pathways.

Broader implications
There are enormous carbon stores sequestered in marinepermafrost in the Arctic, which are of particular concern for release as the warming Arctic Ocean transfers heat to the atmosphere faster than it is transferred from the atmosphere to terrestrial permafrost.Migration from this submerged permafrost reservoir to the ocean has created a vast marine seep field that lies entirely in shallow waters with emissions contributing directly to the atmospheric budget (Shakhova et al., 2014).Widespread ESAS seabed bubble emissions have been documented (Shakhova et al., 2014(Shakhova et al., , 2015) ) demonstrating permafrost integrity failure that makes CH 4 and additional organic carbon available for microbial CH 4 generation.
These observations support the hypothesis that the subsea permafrost is a controlling factor for spatial variability in seabed CH 4 fluxes.The current state of subsea permafrost is key to understanding how CH 4 in ESAS seabed reservoirs escapes to the atmosphere.There is enormous uncertainty in future emissions largely due to the paucity of data.In situ calibrated sonar shows significant promise as a new tool to evaluate seabed fluxes quantitatively over wide areas.

Future directions
In this study, bubble plumes spanning an almost two orders of magnitude flow (0.019 to 1.1 L s −1 ) were studied; however, a key intermediate range (0.045-0.8L s −1 ) was missed.These bubble plumes are in the transition from a nonlinear relationship between σ and flow to a saturation where σ is largely independent of flow.Future experiments should endeavor to follow plumes for more than 15 m; however, currents made this infeasible.Given that seep bubble plumes often escape from nearby vents into plumes that eventually merge and the importance of bubble-bubble acoustic interactions, calibration studies should include multiple bubble plumes from closely located sources.Studies in calmer waters could better elucidate the importance of small bubbles versus large bubbles to overall sonar return.
This study featured the novel use of a numerical bubble plume model to correct for different size evolution between calibration gas bubble plumes and seep bubble plumes.Uncertainty arises from the bubble size distribution, which needs to be measured for the calibration and seep bubble plumes at multiple flow rates.Our approach was a simplified first effort with room for improvement.

Conclusions
In this study, we present a methodology for using an in situ plume calibration approach to derive quantitative sonar methane emissions.We created in situ engineered bubble plumes from a 40 m depth spanning an almost two orders of magnitude flow (0.019 to 1.1 L s −1 ).Nonlinear calibration curves related sonar return to flux for a range of depths and demonstrated significant bubble-bubble acoustic interactions, precluding an inversion approach based on scaling a bubble-sonar cross section with the (unmeasured) size distribution.The weak sonar occurrence probability distribution function was well described by a power law that likely correlated with small-bubble dispersion while strong sonar returns were largely independent of depth, consistent with a focused central core of large bubbles.
The in situ calibration curve was applied to natural seepage from 70 m depth in the Laptev Sea outer shelf where subsea permafrost is predicted to be degraded in modeling studies.A correction then was made for the different volume evolution of the nitrogen calibration plume and the methane seep bubble plume through the use of a numerical bubble plume model.The model was initialized with a typical (assumed) minor bubble plume size distribution and suggested ∼ 5 % correction for the first 5 m depth window.Emissions for three seepage areas of 5.56, 42.73, and 4.88 mmol s −1 were derived from multibeam sonar data with good to reasonable agreement (4-37 %) between single-(biased lower) and multibeam sonar.
The seepage occurrence probability distribution function was bimodal, with weak seepage well described by a power law.This was interpreted as suggesting primarily small minor bubble plumes.The seepage-mapped spatial patterns suggested subsurface geologic control along linear trends.The analysis showed that a probability distribution could provide insights into geologic control.

Figure 1 .
Figure 1.(a) Coal Oil Point (COP) seep field map showing the Shane Seep area of the scoping study.Sonar data from 2005 was adapted from Leifer et al. (2010).(b) Shane Seep multibeam sonar survey map of seep detection (2 m depth window at a seabedfollowing height of 4 m).MBES data was collected in 2009.

Figure 3 .
Figure 3. Locations of oceanographic stations for the RV Viktor Buynitsky cruise, 2012, marked by yellow circles.Polygons of major focus areas are marked as P1 (northern Laptev Sea), P2 (east Lena Delta), and P3 (Dmitry Laptev Strait) shown in insets.Ship tracks accompanied by CTD (conductivity, temperature, and depth) measurements (and geophysical surveys) performed in the P1 are shown as red lines.

Figure 4 .
Figure 4. (a) Salinity and temperature (T ) with respect to depth (z) during engineered bubble plume experiments.(b) Single-beam echosounder sonar return integrated across the plume (σ ) with z for a bubble plume (blue) and for no bubble plume (red); the bubble plume σ is circled.

Figure 5 .
Figure 5. (a) Minor bubble plume size distribution ( ) with respect to the equivalent spherical radius (r) used to initialize the bubble model.(b) Measured temperature (T )-depth (z) profile used in model.

Figure 6 .
Figure 6.Field sonar data from the Coal Oil Point seep field for air bubbles in 22 m deep water.Sonar return counts integrated across the plume (σ ) versus airflow (Q) and height above seabed (h) for four airflows and least-squares linear regression fits to log(σ ) versus h.

Figure 8 .
Figure 8. Sonar return (σ ) versus volumetric flow (Q) calibration curves for the single-beam sonar for (a) all Q and (c) low Q and for the multibeam sonar for (b) all Q and (d) low Q.Fit parameters are shown in TableS2.

Figure 9 .
Figure 9. (a) Depth (z) evolution of the bubble plume size distribution ( ) for a nitrogen minor plume (calibration) from 40 m and (d) for a CH 4 seep plume from 70 m.Seabed normalized volume averaged over depth window (< Q >) of the rising bubble plume for the (b) calibration plume and (e) seep plume.Molar vertical flux for (c) calibration plume, and (f) seep Data keys on panels.

Figure 10 .
Figure 10.Sonar return (σ ) with depth (z) of seep bubble plumes in the Laptev Sea.(a, c, d) Multibeam sonar data, single ping, in each of the seep areas, locations labeled on (b).(b) Single-beam sonar data.The size scales and data keys are shown on the panels.

Figure 11 .
Figure 11.Seep mass flux (Q m ) map for (a) all seep areas and for (b-d) Seep areas 1-3.Data key is shown on panel (c).

Figure 12 .
Figure 12.Seep mass flux (Q m ) occurrence probability distribution function ( (Q m )) normalized to flux bin-width (bin widths are logarithmically-spaced) for (a) all seep areas and for (b-d) Seep areas 1-3 with power law fits.Data key is shown on panel (a).See Table 2 for fits.

Table 2 .
Fit parameters for the seep area flux probability distribution function.