Sonar Estimation of Methane Bubble Flux from Thawing Subsea Permafrost: A Case Study from the Laptev Sea Shelf

: Seeps found o ﬀ shore in the East Siberian Arctic Shelf may mark zones of degrading subsea permafrost and related destabilization of gas hydrates. Sonar surveys provide an e ﬀ ective tool for mapping seabed methane ﬂuxes and monitoring subsea Arctic permafrost seepage. The paper presents an overview of existing approaches to sonar estimation of methane bubble ﬂux from the sea ﬂoor to the water column and a new method for quantifying CH 4 ebullition. In the suggested method, the ﬂux of methane bubbles is estimated from its response to insoniﬁcation using the backscattering cross section. The method has demonstrated its e ﬃ ciency in the case study of single- and multi-beam acoustic surveys of a large seep ﬁeld on the Laptev Sea shelf.


Introduction
Release of previously generated methane (CH 4 ) preserved in natural gas fields, coal beds, and seabed deposits of CH 4 hydrates provides an important feedback in the Arctic climate system. Until recently [1], the Arctic Ocean was not considered as a possible source of CH 4 flux, as the impermeable subsea permafrost was believed [2] to seal marine sediments and prevent methane leakage to the water column and on to the atmosphere from the Arctic seabed which may store significant CH 4 reserves [3,4], including those sequestered in gas hydrates [5]. The permafrost of the East Siberian Arctic Shelf (ESAS), which makes up at least 80% of total subsea permafrost, has trapped the largest hydrocarbon reservoir on the planet [3,6,7]. However, the stability of the sequestered carbon (primarily in CH 4 ) is highly uncertain. Onshore and offshore Arctic permafrost can thaw from the top downward and bottom up. The downward degradation, with the respective expansion of the active layer, produces unfrozen zones (taliks) which can also degrade from below under the effect of geothermal heat flux from underlying unfrozen sediments [8,9]. Thawing from below most often occurs in offshore seabed permafrost [2]. The ongoing warming in the Arctic is especially pronounced over ESAS, where the mean surface air temperature became up to 5 • C higher [10] for the first five years of the 21st century. The permafrost degradation may cause destabilization and dissociation of gas hydrates which can release large amounts of free gas.
The permafrost within ESAS features distinct zones serving as conduits for methane transmission [11] and causing in year-round round methane emission to the atmosphere from the sedimentary reservoir [1,12]. The daily bubble-mediated methane flux from large ESAS seeps can reach hundreds of grams per square meter [12,13]. This makes ESAS a significant modern marine source of CH 4 , contributing to the regional methane budget as much as the terrestrial Arctic ecosystems [11]. The ESAS contribution to atmospheric CH 4 may be even greater, given its area (up to an order of magnitude larger than the Siberian wetlands) and methane emissions throughout winters, when terrestrial ecosystems are dormant [7,12]. Methane emissions in the ESAS are controlled by permafrost degradation, and the future emissions are expected to depend on the coastward dynamics of the subsea permafrost. The amount of CH 4 stored within the shallow ESAS seabed and the apparently widespread thawing of the subsea permafrost suggest that ESAS emissions may increase by 3-5 orders of magnitude [12,14]. The daily bubble-mediated methane flux from large ESAS seeps ( Figure 1) can reach hundreds of grams per square meter [12,13]. Seeps in the Arctic shelf mark the areas of permafrost thawing and related destabilization of gas (methane) hydrates.
top downward and bottom up. The downward degradation, with the respective expansion of the active layer, produces unfrozen zones (taliks) which can also degrade from below under the effect of geothermal heat flux from underlying unfrozen sediments [8,9]. Thawing from below most often occurs in offshore seabed permafrost [2]. The ongoing warming in the Arctic is especially pronounced over ESAS, where the mean surface air temperature became up to 5 °C higher [10] for the first five years of the 21st century. The permafrost degradation may cause destabilization and dissociation of gas hydrates which can release large amounts of free gas.
The permafrost within ESAS features distinct zones serving as conduits for methane transmission [11] and causing in year-round round methane emission to the atmosphere from the sedimentary reservoir [1,12]. The daily bubble-mediated methane flux from large ESAS seeps can reach hundreds of grams per square meter [12,13]. This makes ESAS a significant modern marine source of CH4, contributing to the regional methane budget as much as the terrestrial Arctic ecosystems [11]. The ESAS contribution to atmospheric CH4 may be even greater, given its area (up to an order of magnitude larger than the Siberian wetlands) and methane emissions throughout winters, when terrestrial ecosystems are dormant [7,12]. Methane emissions in the ESAS are controlled by permafrost degradation, and the future emissions are expected to depend on the coastward dynamics of the subsea permafrost. The amount of CH4 stored within the shallow ESAS seabed and the apparently widespread thawing of the subsea permafrost suggest that ESAS emissions may increase by 3-5 orders of magnitude [12,14]. The daily bubble-mediated methane flux from large ESAS seeps ( Figure 1) can reach hundreds of grams per square meter [12,13]. Seeps in the Arctic shelf mark the areas of permafrost thawing and related destabilization of gas (methane) hydrates.
Given the ESAS seepage extent there is a critical need for new effective, rapid, and quantitative monitoring approaches. Sonar surveys provide an effective tool for mapping seabed methane fluxes and monitoring subsea Arctic permafrost seepage driven by permafrost degradation. In this paper we report the results of single-and multi-beam acoustic surveys we undertook for quantifying CH4 ebullition in the ESAS region [12,13].

Materials and Methods
Term gas seep is used in this study to refer to the release of gas in bubbles that rise from the seabed and form stable regions of increased bubble concentration in the water column. Seeps exist in shallow water on the shelf [1,[15][16][17][18] or in deepwater offshore regions [15,[19][20][21]. Each shallow seep occupies a few square km of the sea floor [12,22], while the deepwater seeps are commonly focused within 10 m (point seeps) and dispersed to large distances [15,16,20,21]. The single point seeps, mainly in deep sea, detectable separately by echo sounding are characterized by the flux F, which is the amount of methane carried by rising bubbles from a seep through a horizontal surface per unit time. The densely clustered shallow seeps, which are irresolvable individually, are characterized by the flux Fs from a unit area corresponding to the amount of methane carried by rising bubbles through a horizontal unit area per unit time. Given the ESAS seepage extent there is a critical need for new effective, rapid, and quantitative monitoring approaches. Sonar surveys provide an effective tool for mapping seabed methane fluxes and monitoring subsea Arctic permafrost seepage driven by permafrost degradation. In this paper we report the results of single-and multi-beam acoustic surveys we undertook for quantifying CH 4 ebullition in the ESAS region [12,13].

Materials and Methods
Term gas seep is used in this study to refer to the release of gas in bubbles that rise from the seabed and form stable regions of increased bubble concentration in the water column. Seeps exist in shallow water on the shelf [1,[15][16][17][18] or in deepwater offshore regions [15,[19][20][21]. Each shallow seep occupies a few square km of the sea floor [12,22], while the deepwater seeps are commonly focused within 10 m (point seeps) and dispersed to large distances [15,16,20,21]. The single point seeps, mainly in deep sea, detectable separately by echo sounding are characterized by the flux F, which is the amount of methane carried by rising bubbles from a seep through a horizontal surface per unit time. The densely clustered shallow seeps, which are irresolvable individually, are characterized by the flux Fs from a unit area corresponding to the amount of methane carried by rising bubbles through a horizontal unit area per unit time.

Evaluation of Gas Flux within Water Column Using Sonar: Theoretical Background
Quantifying a methane bubble flux from sonar data is challenging and poorly amenable to laboratory modeling. Currently the problem has been solved by combined theoretical calculations, laboratory testing, and field experiments [12,13,15,[22][23][24][25][26][27][28]. We consider several methods of remote active sonar estimation of bubble methane fluxes based on data from ship-mounted single-and multibeam echosounders operated in a monostatic mode when the acoustic transmitter also acts as a receiver. In this approach, the capacity of a physical object (including bubbles) to reflect and scatter the transmitted sound back to the receiver can be estimated via its backscattering cross section [29], with an area dimension: where I bs is the strength of the scattered signal at the receiver; L is the distance between the transmitter/receiver and the target; I is the strength of the original incident wave. The total (differential) backscattering cross section can be expressed logarithmically as the target strength TS of scatterers in the insonified volume [29,30]: Single bubbles in ebullition zones monitored remotely by echo sounding (sonar) surveys are either resolvable (1) or not (2) in the acoustic data. The flux of methane into water and air due to resolvable single bubbles, which are often quite weak shallow sources [23], can be estimated from their rise velocity and size [15,16,31,32]: where V m is the molar volume of methane at the pressure and temperature of the respective sea depth; S is the insonified area at the given depth; N is the number of bubbles that cross the given surface for the time of observation t; ϑ i is the volume of the i-th bubble [15] or where V b (r) is the volume of a presumably spherical bubble with the radius r; is the averaging operator; P(h) is the pressure at the sea depth h; N b is the number of detected bubbles; t is the observation time; R is the gas constant; T is the Kelvin temperature • K [28]. The insonified area at the given distance L is given by where Ψ D is the integrated beam width [29]: D is the transducer range and dΩ is the solid angle increment. The method is disadvantageous as it requires knowing the bubble size. The common single-beam echosounders cannot pinpoint single bubbles, which leads to large errors. Furthermore, the estimates of highly variable bubble sizes are poorly reliable because bubbles of different sizes can produce the same scattering patterns. Multi-frequency acoustic systems may help overcome this problem, but this solution is beyond the scope of our study.
More accurate backscattering cross-section estimates can be obtained with dual-beam (split-beam) echosounders [33] which use narrow and wide coaxial beams. In this case, calculations include only the scatterers recorded by both beams, and the target is located at the center of the wide beam used for reference. The splitbeam sonars with complex antenna systems and additional phase measurements ensure the best performance [21,24,28,34].
The size of bubbles can be measured optically or inferred from their rise velocity, for instance, using graphical or empirical relationships [12,25,35], including in a special Matlab code [35]. The optical measurement is one of the best method to estimate bubble size, especially using sophisticated camera systems such as Wang et al. (2016) [36]. Nevertheless, these methods demand the chambers in close proximity to the bubbles being studied. Therefore, through it is very difficult, and often impossible, to examine in detail a large area with floating bubbles. However, these estimates of the size of bubbles by the second mentioned method are not very accurate for large bubbles (>1 mm radius) as their size poorly correlates with rise velocity. This velocity can be estimated only if single bubbles are resolvable by the sonar [15], when each bubble appears as an inclined line in the echogram, and the average velocity V is related to its slope as [15]: where ∆h is the depth difference between the ends of the bubble path visible in the image; ∆t is the respective time interval; ϕ is the beam half-width at the level equal to the ratio of signal amplitudes at the path ends to the maximum bubble response. At a high signal/noise ratio, this slope is close to zero beam width, which in common sonars approximately corresponds to 0.7. This equation is valid for the bubble sizes that do not change much for the observation time. The size dependence of the rise velocity [37] can be used to estimate the equivalent spherical radius r i and hence the volume ϑ i of each i-th bubble. The observed maximum height of gas plumes (also called flares) provides another tool for estimating bubble sizes [15,22,27]. It corresponds to the sea depth where the plumes are still seen in echograms as continuous zones of high scattering. If single bubbles in a gas plume are irresolvable but their size, velocity, and shape distributions are known, the methane flux can be inferred from the total back scattering cross section [11,12,15,22,27,38,39], assuming that it equals the sum of cross sections for single bubbles, in the single scattering approximation [40]. In this approximation, the cross section of backscattering from a bubble cloud σ bs is given by [29]: where r is the bubble radius; N 0 is the number of bubbles in the cloud; n(r) is the probability density function of the bubble size distribution; σ bsr is the back scattering cross section of a single bubble of the radius r. The number N 0 depends on the number of bubbles (i) within the h ± l/4 depth interval, where h is the sea depth and l is the pulse length in the case of deepwater gas plumes or single gas sources or (ii) within the segment of integrated beam width (Ψ D ) at a h ± l/4 distance from shallow plumes or plumes rising from a dense cluster of irresolvable sources. Closely spaced bubbles in large gas plumes, where the spacing between bubbles is commensurate with their radius, influence one another [38]. This interaction may reduce markedly the backscattering cross section and lead to underestimation of the respective gas flux values.
The bubble gas flux F (moles per second) is related with the size and shape distribution functions of n(r) and Φ(r), respectively [15,27], as where V(r) is the rise velocity of bubbles.

of 14
Without N 0 in (1) and (2), the gas flux and the backscattering cross section become related as [28,33,41]: where K is a complex distribution function of multiple variables: size, shape, rise velocity, surface, and acoustic properties of bubbles, as well as the sounding frequency, etc. Equation (3) defines the methane flux estimated from backscattering cross sections of deep-sea gas plumes or isolated point sources. In the case of multiple irresolvable sources, the flux from a unit area becomes where S is the effective insonified area. Thus, Equation (3) allows calculating the gas flux using the function K found from the known size, shape, and velocity distributions of bubbles, but these variables are often impossible to constrain. Therefore, the calculations are sometimes simplified by using the average values of the variables, assuming a spherical bubble shape.
For instance, the assumptions in [15] are that all bubbles have the same radius r and rise at the same velocity V; the radius of bubbles exceeds the resonance value; and kr << 1, where k is the wavenumber. With these assumptions, the methane flux (mole/s) is given by where V M is the molar volume of methane at the temperature and pressure within the sounding depth, and l is the pulse length. Otherwise, some variables are used as average values and others as approximating relationships, such as the exponential approximation of the bubble size distribution n(r) = e −∝r [28] that decreases stepwise to zero at r ≤ r c (r c is the smallest detectable size of bubbles; r ≥ r c and n(r) = 0 at r ≤ r c ). According to laboratory and field data, r c is from 0.5 to 1.0 mm and α = 1.1 mm −1 .
Note that the radii of the rising bubbles estimated as above are tentative. The size, shape, and velocity distributions of bubbles require more rigorous constraints, for example, by optical methods that allow direct measurements of these parameters for each bubble [12,22,25,35].
In all above solutions for the methane bubble flux, the backscattering cross section of a single bubble (σ bsr ) was found as where r res is the resonance radius of the bubble at the sonar operation frequency and δ is its damping factor. This is the most widespread approach in the sonar gas flux estimation practice [12,14,22,23,28]. Note that Equation (4) is valid at kr 1, which imposes limitations on the sounding frequency. Strictly speaking, the 38 kHz frequency often used for methane flux estimation [21,28] is too high, because kr = 1 already at r = 6 mm. However, the validity of (5) can be extended into the kr > 1 domain to a satisfactory accuracy if the real nonspherical shape of bubbles is taken into account [27]. On the other hand, according to (5), the backscattering cross section of bubbles smaller than the resonance radius, decreases with decreasing bubble size at r 6 , i.e., the contribution of these bubbles to the backscattered signal is vanishing, which is especially essential for deep-sea seeps. The resonance radius of bubbles is given by [27]: where f 0 is the operation frequency; γ is the adiabatic exponent for gas in the bubble; P 0 is the atmospheric pressure; ρ w is the water density; g is the acceleration due to gravity; h is the sea depth. The resonance radii estimated with Equation (6) are 0.6 and 0.4 mm at a depth of 50 m (reached by the common 12 and 18 kHz sonar surveys), and 2.6 and 1.6 mm at 1 km. For shallow seeps, the frequency 18 kHz satisfies both conditions kr 1 and r c ≥ r res , as one can infer from the observed bubble sizes [12,25], though other frequencies, including hundreds of kHz, are applicable as well. The high-frequency systems are advantageous by being small, mobile, and energy-saving and thus convenient for research on various ships, even those that lack the mounted specialized equipment.
Choosing the sonar operation frequency for the responses of deep-sea seeps is more difficult as the range fitting both r c ≥ r res and kr 1 decreases with depth, which may lead to large errors in the gas flux estimates.
Note that shallow seeps may pose additional problems due to strong resonance effects, as the attenuation Q = 1/δ of bubbles increases at decreasing pressure. The damping factor cannot exceed 7 at sea depths about 1000 m, where it is limited by the theoretical value associated with re-radiation, but it can reach a few tens at shallower <100 m depths [42]. Thus, the resonance effects increase markedly in the shallow sea, and the backscattering cross section of a single resonant bubble can be larger than that of several non-resonant ones.
Note that the simplifications used to estimate the bubble methane flux from backscattering cross section of a single bubble as in (5), with simplified approximations or average values instead of real size and rate distributions, are not always valid and may cause significant uncertainty. Furthermore, it is often impossible to choose the frequency fitting both kr 1 and r c ≥ r res conditions, especially in the case of deep-sea seeps. As for the responses of shallow seeps, a large error may result from poor accuracy of the Q factor of bubbles and approximated attenuation of acoustic waves, as well as from neglected shape distribution of bubbles. Therefore, it is important to obtain rapid estimates of K in (3) in the conditions most similar to the real field surveys. Since this coefficient is a complex function of multiple variables, which are mostly unknown, the methods for obtaining its empirical rather than calculated values [12,13,22,34] are of special interest. So in work [34] on the basis of experiments with a horizontally looking single beam transducer (40 and 300 kHz) directed towards artificially produced bubbles, it was shown that an acoustic system can be calibrated in such a way that gas flux rates of bubble-size spectra, as observed at natural seeps, can be directly related to the echo level of a known, acoustically insonified volume.

Evaluation of Gas Flux within Water Column Using Sonar: Field Calibration Method
We estimated the methane bubble flux from the sea floor into the water column by a method similar to that in [34]. The main advantage of this approach is that the gas flux estimation does not require the knowledge of size, shape, and velocity distributions of bubbles and the specifications of the echosounder. Our method [12] differs from that suggested in [34] in two main aspects: the calibration in the field rather than in laboratory and vertical rather than horizontal insonification. The calibration was performed in typical local field conditions, while the vertical transmission and reception of signals was chosen since the methane flux was measured bottomward by ship-based echosounders. The measurements were run in the monostatic mode, as in [34], when the sonic waves are emitted and received by the same electroacoustic transducer. Since the transducer was submerged a few meters below the water surface, the distance L to the scattering volume was slightly less than its depth h. At the depth h, the echo is produced by rising bubbles within the effective scattering volume ∆V, which for h cτ/2 is given by [29]: where c is the sound speed; τ is the pulse duration; Ψ D is the integrated beam width. The target strength within the effective scattering volume ∆V is described by the standard sonar equation [30]: where TS is the target strength of bubbles within the effective scattering volume ∆V; RL is the received level, in dB/mPa; SL is the source level in the transmitter, in dB/mPa at 1 m; TL is the transmission loss due to spherical divergence and attenuation: where L 1 is the distance between the acoustic projector and the scattering bubbles (1 m) and α is the attenuation coefficient of the seawater, in dB/m. The received level is found as where U 2 n is the normalized amplitude of the received signal for each ping n; K 1 is the constant depending on the parameters of the echosounder.
Commonly the modern sonars are designed to compensate the transmission losses. Thus, without TL and with constant SL, the TS equation becomes where K 2 is the calibration coefficient; U 2 0 is the normalized amplitude of each ping n at compensated transmission losses. Therefore, the methane flux from isolated seeps is The flux from unit area in shallow seeps, with acoustically irresolvable bubbles, is where S is the effective scattering area. As follows from (8), the calculations can be limited to simple calibration measurements of gas flux dependence of the backscatter strength for deep-sea seeps, but additional knowledge of angle parameters may be required to estimate the effective scattering area S by (9) in the case of irresolvable shallow seeps. This area is found from the distance between the projector and the insonified volume and from the integrated beam width Ψ D . Unfortunately, the manufacturer does not provide Ψ D values in the specifications. For a piston transducer, it can be estimated as [29] Ψ D = 5.48 (kr) 2 where k is the wavenumber; r is the transmitter radius. Analysis of multibeam records is worth special interest as they include backscattering strength in 3D, depth and period of emitted signals, as well as beam number of the multibeam system. The fragment of a multibeam sonar cross-section of a seep in Figure 2a shows distortions associated with propagation of a sonic wave that was transmitted and received at different angles to the vertical in the water column. In other words, a 90 • incident signal will reach the 100 m depth before that emitted at 60 • . The shift of the signal in depth (L) and insonified area (D) can be calculated from the known sound speed (c), sampling rate ( f qua ), beam number (n) and dip (α), and the distance (n) traveled by the signal, using Equations (10) and (11).
where c is the sound speed; f qua is the sampling rate; n is the beam number; α is the beam dip; i is the distance traveled by the signal.
Geosciences 2020, 10, x FOR PEER REVIEW 8 of 14 The variables and , as well as the beam number , should be added to Equations (4) and (9) for the bubble gas flux estimation.

Calibration of the Sonar
The survey campaign at a test site in the East Siberian Arctic shelf included calibration work in accordance with the project objectives: (i) standard procedure using a tungsten carbid target 38.1 mm in diameter and (ii) special calibration using an artificial sounder-derived gas plume. The data were processed in MatLab, with estimation of mean values and standard deviation. The match between different data populations was checked with the Student criterion (p < 0.05).
The standard calibration consisted in estimating the maximum square amplitude of the received signal that was backscattered by a target submerged to 4.5 m below the water surface in the middle of the insonified area, at user-specified sounding parameters (pulse frequency, amplitude, duration, etc.). The normalized amplitude represents each signal with compensated transmission losses, proceeding from the real distribution of temperature, salinity, and sound speed in the given place of the water column. The conductivity, temperature, and depth data were acquired with a CastAway-CTD probe within 100 m depths, at a sampling rate of 5 Hz, to an accuracy of ±0.1 PSU for salinity and ±0.05 °C for temperature. The strength of the calibration target was estimated as in [43][44][45] from temperature and salinity at the target depth (Figure 3a). The variables L and D, as well as the beam number n, should be added to Equations (4) and (9) for the bubble gas flux estimation.

Calibration of the Sonar
The survey campaign at a test site in the East Siberian Arctic shelf included calibration work in accordance with the project objectives: (i) standard procedure using a tungsten carbid target 38.1 mm in diameter and (ii) special calibration using an artificial sounder-derived gas plume. The data were processed in MatLab, with estimation of mean values and standard deviation. The match between different data populations was checked with the Student criterion (p < 0.05).
The standard calibration consisted in estimating the maximum square amplitude of the received signal U 2 0 that was backscattered by a target submerged to 4.5 m below the water surface in the middle of the insonified area, at user-specified sounding parameters (pulse frequency, amplitude, duration, etc.). The normalized amplitude U 2 0 represents each signal with compensated transmission losses, proceeding from the real distribution of temperature, salinity, and sound speed in the given place of the water column. The conductivity, temperature, and depth data were acquired with a CastAway-CTD probe within 100 m depths, at a sampling rate of 5 Hz, to an accuracy of ±0.1 PSU for salinity and ±0.05 • C for temperature. The strength TS of the calibration target was estimated as in [43][44][45] from temperature and salinity at the target depth (Figure 3a).
According to the frequency dependence of TS obtained by calibration of the Simrad EK15 echosounder (Figure 3b), the average TS is −39.18 dB/m 2 , which corresponds to the σ bs0 = 1.21 cm 2 backscattering cross section of the calibration sphere.
The calibration provided constraints on the correlation K 0 between the target scattering cross section σ bs0 and the normalized square amplitude of the calibration signal: Geosciences 2020, 10, 411 9 of 14 Special calibration of the Simrad EK15 echosounder was performed using a gas plume produced by a bubble generator [12,13,22] consisting of several sequentially connected units: a gas bomb with nitrogen, a gas supply system with gas flux control, and a submerged nozzle with an outlet diameter of 3 mm; the flux rate was set in a range of 0.02 L/s to 1.27 L/s. processed in MatLab, with estimation of mean values and standard deviation. The match between different data populations was checked with the Student criterion (p < 0.05).
The standard calibration consisted in estimating the maximum square amplitude of the received signal that was backscattered by a target submerged to 4.5 m below the water surface in the middle of the insonified area, at user-specified sounding parameters (pulse frequency, amplitude, duration, etc.). The normalized amplitude represents each signal with compensated transmission losses, proceeding from the real distribution of temperature, salinity, and sound speed in the given place of the water column. The conductivity, temperature, and depth data were acquired with a CastAway-CTD probe within 100 m depths, at a sampling rate of 5 Hz, to an accuracy of ±0.1 PSU for salinity and ±0.05 °C for temperature. The strength of the calibration target was estimated as in [43][44][45] from temperature and salinity at the target depth (Figure 3a).  Calibration provided an empirical relationship between the backscattering strength at the depth 39 m and the gas flux from a nozzle at the depth 40 m in the middle of the beam width. Since any gas other than methane can be used for calibration [13], we chose nitrogen, which is a safe inert gas. The diameter of bubbles estimated by direct optical measurements ranged from 3 to 12 mm and corresponded to the characteristic size of natural bubbles in seep sites [22,25,31].

Relationship between Gas Flux and Measured Response
The difference between the calibrations when using methane and nitrogen can be neglected, since for the densities of these gases of the used depths (39 m) are much lower than the density of water. For bubbles, the difference at different gas content is observed in the region of resonance frequencies.
Since the operating frequencies are much higher than the resonant frequencies, this difference can be neglected. The difference in diffusion of these two gases can also be neglected, since measurements are taken near seabed. In addition, these depths are much higher than the zone of stability of methane gas hydrate (about 300 m), so the formation of a gas hydrate crust on the surface of the bubbles does not occur. The applied 200 kHz frequency was much higher than the resonance frequency of such bubbles and thus satisfied the calibration requirement. The calibration sonar cross-section of Figure 4a was obtained at a flow rate of 1.27 L/s and resolved well the sea floor and the zone of ebullition-related scattering in the water. The signal attenuated at depths shallower than~15 m due to a strong near-surface current which dragged the rising bubbles outside the beam width.
The calibration curve (Figure 4b) represents the gas flux dependence of normalized square voltage at the echosounder output, with each point being a result of averaging over 420 pings in 7 min sampling intervals. The experimental points fit well the straight line through the origin of coordinates which records linear correlation between the signal (U 2 0 ) and the flux (F) in the used range from 1 to 57 mmol/s and thus confirms the validity of Equations (8) and (9) for a large range of gas fluxes. Therefore, the empirical coefficient K obtained in field conditions can make basis for rapid estimation of bubble methane fluxes without data on bubbles in gas plumes and echosounder specifications. This calibration method is applicable to multibeam sounders as well [9]. stability of methane gas hydrate (about 300 m), so the formation of a gas hydrate crust on the surface of the bubbles does not occur. The applied 200 kHz frequency was much higher than the resonance frequency of such bubbles and thus satisfied the calibration requirement. The calibration sonar crosssection of Figure 4a was obtained at a flow rate of 1.27 L/s and resolved well the sea floor and the zone of ebullition-related scattering in the water. The signal attenuated at depths shallower than ~15 m due to a strong near-surface current which dragged the rising bubbles outside the beam width.  Note that in our case, the calibration was carried out for a source at a depth of 40 m. These data were extrapolated to other depths by introducing the values of the molar volume V M and sonicated area into the Formula (9), recalculated for a given depth and temperature at the bottom.

Results and Discussion
Sonar estimation of bubble-mediated methane fluxes requires an operation frequency f exceeding the resonance frequency of rising bubbles f r but low enough to reduce the damping of transmitted and echoed signals. We used a 200 kHz Simrad EK15 singlebeam echosounder with its frequency above the resonance values of bubbles at all selected sea depths from 0 to 100 m. This study covered area of the ESAS: the ice-free area of the Laptev Sea (between 76.5-77.5 • N and 121-132 • E, water depth between 50 and 165 m ( Figure 5).
Geosciences 2020, 10, x FOR PEER REVIEW 10 of 14 range from 1 to 57 mmol/s and thus confirms the validity of Equations (8) and (9) for a large range of gas fluxes. Therefore, the empirical coefficient obtained in field conditions can make basis for rapid estimation of bubble methane fluxes without data on bubbles in gas plumes and echosounder specifications. This calibration method is applicable to multibeam sounders as well [9].
Note that in our case, the calibration was carried out for a source at a depth of 40 m. These data were extrapolated to other depths by introducing the values of the molar volume and sonicated area into the Formula (9), recalculated for a given depth and temperature at the bottom.

Results and Discussion
Sonar estimation of bubble-mediated methane fluxes requires an operation frequency exceeding the resonance frequency of rising bubbles but low enough to reduce the damping of transmitted and echoed signals. We used a 200 kHz Simrad EK15 singlebeam echosounder with its frequency above the resonance values of bubbles at all selected sea depths from 0 to 100 m. This study covered area of the ESAS: the ice-free area of the Laptev Sea (between 76.5-77.5° N and 121-132° E, water depth between 50 and 165 m ( Figure 5). The methane flux quantified from measured backscattering cross sections of bubbles in a large methane seep from the Laptev shelf [12] was compared with the results of calibration against an engineered gas plume. The respective ebullition zone is clearly pronounced in the example sonar cross-section of Figure 6a.
Echo sounding surveys of 49.7 km over a test area of 6 km 2 2018 in a region of a large methane seep discovered in 2011 [12,46] revealed 46 independent vents from 15 to 336 m in diameter, including 15 seeps larger than 100 m (Figure 6b). The seeps were found out to cover 1.2 km 2 and produce a gas The methane flux quantified from measured backscattering cross sections of bubbles in a large methane seep from the Laptev shelf [12] was compared with the results of calibration against an engineered gas plume. The respective ebullition zone is clearly pronounced in the example sonar cross-section of Figure 6a. We also performed observations in the southernmost part of the Laptev Sea, in Ivashkina Lagoon, which has been progressively inundated during the last ~100-500 years, replacing a former thermokarst lake. Bubble release occurred from narrow, steep depressions aligned parallel to the lagoon's northern edge. Backscattering cross-sections of the bubbles emitted from 17 seeps observed in Ivashkina Lagoon were recorded for 36 h using portable single-beam sonar, which was calibrated in situ during the same campaign. In Ivashkina Lagoon, CH4 fluxes observed in October 2013 ranged from 5 to 24 g m −2 d −1 [12].
Thus, the proposed method, based on the determination of the empirical coefficient K in field conditions, makes it possible to carry out a quick estimation of bubble methane fluxes without data on bubbles (distributions of size, shape and rise velocity) in gas plumes and echo sounder calibrations. This method is applicable to multibeam echo sounders as well.

Conclusions
Among several other reviewed approaches, the new method for sonar quantifying of methane ebullition in the water column based on calculations from the cross section of backscattering from rising СН4 bubbles has been tested in a large seep area on the Laptev Sea shelf has demonstrated its efficiency. The gas flux estimates by the new method are comparable to those obtained by calibration against an engineered gas plume: 5.7 ± 0.6 and 6.8 ± 0.8 kg/s, respectively. Thus, both methods are applicable for rapid remote estimation of СН4 fluxes from seeps, taking into account that the backscattering cross section may be about 20% lower than the true values.  Echo sounding surveys of 49.7 km over a test area of 6 km 2 2018 in a region of a large methane seep discovered in 2011 [12,46] revealed 46 independent vents from 15 to 336 m in diameter, including 15 seeps larger than 100 m (Figure 6b). The seeps were found out to cover 1.2 km 2 and produce a gas flux of 5.7 ± 0.6 kg/s, as estimated from measured backscattering cross section and by calculations with Equation (4), using measured radius and rise velocity of bubbles; the >100 m seeps provide the greatest contribution to the total flux (Figure 6c). The estimates by Equation (9) using data calibrated against an artificial gas plume gave 6.8 ± 0.8 kg/s. The two gas flux values differ to 19%, at p < 0.05.
We also performed observations in the southernmost part of the Laptev Sea, in Ivashkina Lagoon, which has been progressively inundated during the last~100-500 years, replacing a former thermokarst lake. Bubble release occurred from narrow, steep depressions aligned parallel to the lagoon's northern edge. Backscattering cross-sections of the bubbles emitted from 17 seeps observed in Ivashkina Lagoon were recorded for 36 h using portable single-beam sonar, which was calibrated in situ during the same campaign. In Ivashkina Lagoon, CH 4 fluxes observed in October 2013 ranged from 5 to 24 g m −2 d −1 [12].
Thus, the proposed method, based on the determination of the empirical coefficient K in field conditions, makes it possible to carry out a quick estimation of bubble methane fluxes without data on bubbles (distributions of size, shape and rise velocity) in gas plumes and echo sounder calibrations. This method is applicable to multibeam echo sounders as well.

Conclusions
Among several other reviewed approaches, the new method for sonar quantifying of methane ebullition in the water column based on calculations from the cross section of backscattering from rising CH 4 bubbles has been tested in a large seep area on the Laptev Sea shelf has demonstrated its efficiency. The gas flux estimates by the new method are comparable to those obtained by calibration against an engineered gas plume: 5.7 ± 0.6 and 6.8 ± 0.8 kg/s, respectively. Thus, both methods are applicable for rapid remote estimation of CH 4 fluxes from seeps, taking into account that the backscattering cross section may be about 20% lower than the true values.