What the bubble knows: Lake methane dynamics revealed by sediment gas bubble composition

Atmospheric methane (CH4) concentrations have more than doubled in the past ~ 250 yr, although the sources of this potent greenhouse gas remain poorly constrained. Freshwaters contribute ~ 20% of natural CH4 emissions, about half attributed to ebullition. Estimates remain uncertain as ebullition is stochastic, making measurements difficult, time consuming, and costly with current methods (e.g., floating chambers, funnel gas traps, and hydroacoustics). We present a novel approach to quantify basin‐wide hypolimnetic CH4 fluxes at the sediment level based on measurements of bubble gas content and modeling of dissolved pore‐water gases. We show that the relative ebullition flux pathway can be resolved by knowledge of only bubble gas content. As sediment CH4 production, diffusion, and ebullition are interrelated, the addition of a second observation allows closing the entire sediment CH4 balance. Such measurements could include bubble formation depth, sediment diffusive fluxes, ebullition, sediment CH4 production, or the hypolimnetic CH4 mass balance. The measurement of bubble gas content is particularly useful for identifying local ebullitive hotspots and integrating spatial heterogeneity of CH4 fluxes. Our results further revealed the crucial effect of water column depth, production rates, and hypolimnetic dissolved CH4 concentrations on sediment CH4 dynamics. Although we apply the model to cohesive sediments in an anoxic hypolimnion, the model can be applied to shallow, oxic settings by altering the CH4 production rate curve to account for oxidation. Utilizing our approach will provide a deeper understanding of in‐lake CH4 budgets, and thus improve CH4 emission estimates from inland freshwaters at the regional and global scales.

Atmospheric methane (CH 4 ) concentrations have more than doubled in the last 250 yr (Wuebbles and Hayhoe 2002). With a global warming potential 34 times higher than CO 2 over a 100 yr period (Myhre et al. 2013), CH 4 has come to renewed attention as atmospheric concentrations continue to rise after plateauing during the last decade (Dlugokencky et al. 2011;Nisbet et al. 2014). There is currently considerable uncertainty in the global atmospheric CH 4 budgets and the reasons behind the observed atmospheric CH 4 trends, highlighting the need to better understand the various sources and sinks (Kirschke et al. 2013;Turner et al. 2017;Worden et al. 2017). Lakes as natural sources of microbially produced CH 4 contribute about 20% of the natural emissions to the global budget (Bastviken et al. 2011). Approximately, half of these emissions are attributed to CH 4 ebullition; however, these estimates remain highly uncertain (Bastviken et al. 2011) as ebullitive emissions are notoriously difficult, time consuming, and costly to accurately quantify (Wik et al. 2016).
In freshwater sediments, significant quantities of biogenic CH 4 are naturally produced by microorganisms (methanogens) in the terminal process of anaerobic degradation of organic matter (Conrad 2005). Two main pathways are distinguished: acetoclastic methanogenesis utilizing acetate and hydrogenotrophic methanogenesis, which is based on the reduction of CO 2 with hydrogen (Conrad 2005). CH 4 production rates in sediments have been found to be dependent on the quality of organic matter, O 2 exposure, and temperature (Schulz and Conrad 1995;Sobek et al. 2012;Yvon-Durocher et al. 2014). Production rates generally decrease rapidly and near exponentially with sediment depth (Zepp Falz et al.  Popp et al. 2000;Wilkinson et al. 2015). If the net CH 4 production is higher than the diffusive transport from the sediment to the water, sediment pore-water dissolved gas concentrations will become oversaturated and bubbles are formed and may be released (Boudreau 2012;Schmid et al. 2017).
For pore-water bubble formation, it is necessary to consider the total dissolved gas pressure, which is the sum of the dissolved partial pressures of all sparingly soluble sediment gases (mainly CH 4 and nitrogen, N 2 ). Bubbles form when the total dissolved gas pressure becomes equal to or exceeds the local pressure (hydrostatic + atmospheric; Miyake 1951; D'Aoust 2007; see Fig. 1). Even though CO 2 is biologically produced in the sediment (Conrad 2005), total dissolved gas pressure is mainly driven by CH 4 as CH 4 is~29 times less soluble than CO 2 (at 5 C) (Sander 2015). Accordingly, sediment bubbles contain little CO 2 and consist mainly of CH 4 and N 2 (Casper et al. 2000;Walter et al. 2008). The mechanics of formation, growth, and release of sediment bubbles (ebullition) depend on physical properties of the sediments (Boudreau 2012;Katsman et al. 2013;Liu et al. 2016Liu et al. , 2018Scandella et al. 2016) and changes in local pressure and bottom shear stress (Casper et al. 2000;Joyce and Jewell 2003).
Ebullitive release of sediment gas is highly variable in space and time (DelSontro et al. 2015). Consequently, considerable efforts are required to obtain accurate estimates for whole-lake emissions (Natchimuthu et al. 2016). CH 4 ebullition is typically measured using floating chambers (Bastviken et al. 2004), with hydroacoustic methods (Ostrovsky et al. 2008) or with gas trap funnels (DelSontro et al. 2010;Delwiche et al. 2015;Wilkinson et al. 2015). These conventional approaches have several limitations. First, due to the stochastic nature of ebullition, a large number of measurement sites are required to obtain accurate estimates representative for a whole lake. With gas trap funnels, for example, Wik et al. (2016) recommended 39 measurements of ebullition at 11 or more locations, which then provides an unbiased estimation of ebullition. Second, both gas trap funnel and hydroacoustic methods rely on the knowledge of CH 4 fraction in bubble gas to determine CH 4 ebullitive rates. Even though it is possible to analyze the CH 4 content in trapped bubble gas, measurements may be biased due to gas exchange with the water column and potential CH 4 oxidation. Finally, no method is available that integrates all components needed for a mechanistic understanding of sediment CH 4 mass balance (consisting of production, diffusion, and ebullition). To overcome these challenges, we propose a novel approach based on modeling of dissolved pore-water gases combined with knowledge of CH 4 bubble content to improve basin-wide CH 4 dynamics, including ebullitive emission estimates.
Several authors have demonstrated that CH 4 ebullition deprives pore water of other gases, some of which would otherwise reflect atmospheric saturation concentrations (e.g., N 2 and noble gases ;Reeburgh 1969;Martens and Berner 1977;Kipphut and Martens 1982;Chanton et al. 1989;Brennwald et al. 2005). Reeburgh (1969) investigated dissolved gases in sediments from active ebullition sites and found profiles of N 2 and argon (Ar) to be depleted with respect to concentrations at the sediment-water interface (SWI). The author attributed the phenomenon to bubbles transporting both CH 4 and other pore-water dissolved gases (i.e., N 2 and noble gases) out of the sediment (Reeburgh 1969;Brennwald et al. 2005). This effect, which Reeburgh (1969) termed "stripping," is primarily dependent on ebullition rate and water column depth. Both Martens and Berner (1977) and Kipphut and Martens (1982) estimated ebullition rates by inverse modeling of pore-water profiles of dissolved N 2 and Ar. In agreement with the stripping hypothesis, Chanton et al. (1989) observed that at higher ebullition rates, the bubble gas leaving the sediment contained less N 2 and more CH 4 . Chanton et al. (1989) suggested that this information could be useful for estimating ebullitive fluxes.
More recently, Bazhin (2003Bazhin ( , 2010) developed a steady state theory for dissolved pore-water gases in sediments with active CH 4 ebullition, which has been verified in the laboratory by Kusmin et al. (2006). Bazhin (2003Bazhin ( , 2010 introduced an upper diffusive layer below the SWI where bubbles do not form since the total dissolved gas pressure does not reach oversaturation. The model is driven entirely by sediment CH 4 production, which can be parameterized by an exponential decay. From the exponential fit of sediment CH 4 production Fig. 1. Schematic of CH 4 processes for pore-water mass balance model with the ebullition model based on Bazhin (2003Bazhin ( , 2010. (a) Sediment CH 4 production, (b) CH 4 diffusion flux F CH4,Diff at SWI, (c) total ebullition flux F Ebul,tot at the SWI, and (d) composition of sediment bubble gas. P i is the dissolved gas partial pressure of a given gas. with depth, the model calculates the thickness of the diffusive layer, pore-water profiles of dissolved CH 4 and N 2 , CH 4 diffusive and ebullitive fluxes, and finally the bubble gas composition (i.e., content of CH 4 and N 2 ; see Fig. 1). The model of Bazhin (2003Bazhin ( , 2010 has two degrees of freedom. Therefore, the state of the system, including the complete sediment CH 4 mass balance, can be determined from two given independent observations (discussed below). CH 4 produced in the sediment (Fig. 1) can exit via diffusion and ebullition. In a steady state situation, the mass balance at the sediment level is given by net production = diffusion + ebullition. In this study, we present an inverse modeling approach that utilizes this relationship combined with sediment bubble gas composition (Fig. 1d) to estimate CH 4 pathways in the small, eutrophic kettle lake, Lake Soppen (aka Soppensee).
Specially, we apply a sediment pore-water model and show that: • With reasonable assumptions, the fraction of the produced CH 4 leaving the sediment either as ebullition or diffusion can be predicted with only knowledge of the bubble CH 4 content and the water depth where the bubble was collected.
• With the addition of any additional independent observation (below), the complete sediment CH 4 mass balance can be resolved.
The additional observation can be any of the following: (1) point ebullition estimates, (2) depth of sediment-bubble formation, (3) a well-resolved sediment CH 4 production profile, or (4) CH 4 diffusive flux. Given the same measurement effort, traditional methods (i.e., funnels, hydroacoustics, surface chambers, micrometeorological, open path optical, and infrared camera systems) combined with our approach provide significantly improved estimates of whole-lake CH 4 dynamics.

Study site
Lake Soppen (aka Soppensee) is a small, glacially formed kettle lake in Buttisholz (Canton Lucerne) in the Swiss plateau at an elevation of 596 m above sea level. It has a surface area of 0.227 km 2 , a maximal and mean depth of 26 and 12.3 m ( Fig. 2; Lotter 1989). Given the small catchment area of 1.6 km 2 and a mean annual outflow of 0.03 m 3 s −1 through a small stream in the north, Fischer (1996) estimated the maximum residence time of the water to be 3.1 yr. The lake is eutrophic and lies in an area of intense agriculture. The sediments in the top 6 m are authigenic, fine-grained cohesive, and rich in organic carbon, whereas below 6 m, glacially deposited clays are found (Lotter 1991;Fischer 1996). From the observed depletion of dissolved noble gases in the sediment pore water (down to 7 m), Brennwald et al. (2005) estimated that the lake sediments have been releasing CH 4 bubbles for the past centuries.
The lake was monitored from April 2016 to December 2017 for CH 4 and water column properties, including temperature, O 2 , and CH 4 in the water column (see Supporting Information section S1 and example profile in Supporting Information Fig. S1). In summer, water temperature reaches 24-26 C in the surface mixed layer, whereas temperatures remain close to 5 C in the bottom waters below 10 m water depth. During the stratified periods, the lake water becomes entirely anoxic below~8 m, whereas dissolved CH 4 concentrations increase rapidly below 10 m, becoming as high as 1.3 mol m −3 at the bottom of the lake in autumn 2017.

Sediment CH 4 flux estimates and measurements
We quantify the three components of the steady state sediment CH 4 mass balance: (1) sediment CH 4 production, (2) diffusive CH 4 flux at the SWI, and (3) ebullitive CH 4 flux at the SWI. Sediment CH 4 production was measured from the increase of CH 4 in incubation vials filled with sediment subsamples from four cores taken at water depths 8, 20, and 26 m on the 20 July 2016 and 26 m water depth on the 01 March 2016 (see Supporting Information Figs. S2-S5 and Tables S1-S3). As in situ and incubation temperatures were slightly different, we applied a temperature correction based on the data of Nozhevnikova et al. (2007) (see Supporting Information   S6). The sediment incubation methods followed the procedure by Wilkinson et al. (2015) and are detailed in the Supporting Information section S2.
Diffusive CH 4 fluxes at the SWI were calculated from porewater measurements of cores sampled from various depths in 2016 and 2017. The measurements followed the methods of Donis et al. (2017) and are described in details in the Supporting Information section S3. Briefly, a cut syringe method was applied and head space gas was analyzed on a Cavity Ringdown Spectrometer (Model Picarro G2201-i), which provided both pore-water CH 4 concentration and δ 13 C CH4 values (the Picarro instrument came with factory calibration and a comparison with a brand-new instrument [same model] showed excellent agreement for the δ 13 C CH4 values). Diffusive CH 4 fluxes at the sediment-water interface were then calculated according to Fick's first law (see Supporting Information Fig. S7 and Table S4).
Total ebullitive and CH 4 ebullitive fluxes at the SWI were measured with inverted funnel gas traps at three locations in the lake with water depths 10, 15, and 21 m (see Fig. 2). The funnels were installed above the sediment at 4 m below the water surface and captured gas bubbles released from the sediment. With the bubble model of McGinnis et al. (2006), we subsequently corrected fluxes for gas exchange during the rise of bubbles from the SWI to the funnel gas trap and thus obtained rates of total gas ebullition and CH 4 ebullition at the SWI. We show ebullition rates from the period July to December 2017, during which data are available for all the three funnel sites (see Supporting Information section S4 and Table S5).

Sediment bubble gas sampling, composition, and origin
We collected bubble gas samples from the sediments withiñ 5 m of the funnel sites and along a transect (see Fig. 2). The sampler consisted of a funnel with a weight attached underneath (i.e., at the wide part) and a glass crimp vial (50 or 120 mL) attached to the top (see Supporting Information section S5 and sketch of sampling device in Supporting Information Fig. S8). To sample the bubble gas, the device was filled with lake water and lowered to the sediments. The sediment was impacted~15 times with the weight to completely fill the vial with gas. When the sampling device was retrieved to lake surface, the vial was capped (polytetrafluoroethylene-coated butyl rubber) and crimped under water. Sample vials were kept submerged in lake water until measurement within 1 d.
The composition of the gas samples (CH 4 , CO 2 , N 2 , and O 2 ) was determined with a field portable mass spectrometer (see Brennwald et al. [2016] and Supporting Information section S6 and Table S6). We measured the stable carbon isotopic signature (δ 13 C CH4 ) of bubble gas samples on the same Picarro stable isotope analyzer used for pore-water analysis. Samples were diluted to concentrations in the range of the stable isotope analyzer (< 500 ppm) by injecting 50 μL of bubble gas using a 0.5-mL gas-tight glass syringe (SGE Analytical Science) into a 120 mL sealed glass vial containing artificial air (no CH 4 or CO 2 , Carbagas). After injection of the sample, an additional 60 mL of artificial gas was injected into the 120 mL glass vial, with another 60 mL syringe inserted to relieve the pressure. The two inserted 60-mL all-glass syringes (Poulten Graf Fortuna Optima) were then used to mix the gas in the 120 mL vial by alternatively pressing and releasing the plungers. This was done five times to ensure adequate mixing.
The sediment depth of bubble origin was estimated based on the assumption that the δ 13 C CH4 of the dissolved porewater CH 4 should match the δ 13 C CH4 of the bubble gas CH 4 . For this, we calculated the average δ 13 C CH4 pore-water profile obtained between 8-13 m (n = 5) and 18-26 m water depth (n = 5) and took their average for the depth range between 13 and 18 m water depth (no pore-water profiles were measured in that depth range). We subsequently utilized these δ 13 C CH4 profiles to infer the bubble origin depth from δ 13 C CH4 of the sediment bubble gas as shown for an example in Supporting Information Fig. S9 and section S7.

Model introduction
Our modeling approach largely follows the theory of dissolved gases in sediment pore water as described in Bazhin (2003Bazhin ( , 2010. The model describes the concentrations of dissolved gases CH 4 and N 2 (and potentially other gases) in the sediment pore water (see overview in Fig. 3). Except for CH 4 , which is produced, the pore-water dissolved gas concentration profiles are assumed to be exclusively driven by molecular diffusion and transport of stripped pore-water gases out of the sediment by ebullition (i.e., bubble stripping). The profiles of dissolved gas concentration are assumed to be in steady state and that no further gas transport occurs across the bubble surface once the bubble moves from its depth of origin. As in Bazhin (2003Bazhin ( , 2010, the upper sediment diffusive layer, where theoretically no ebullition occurs, is defined where the total dissolved gas pressure is lower than the local (atmospheric + hydrostatic) pressure (Fig. 3).
The model development is described in the next section, with the key variables listed in Table 1 and the principle assumptions summarized below. The principle assumptions should be considered when applying the model; we revisit the major assumptions in the discussion.
• Conceptualization with upper nonebullitive layer (diffusion layer) and lower ebullition layer (see Bazhin 2010 and assumptions therein). • Dissolved pore-water concentration profiles are in steady state and the average seasonal effect of intermittent ebullition on pore-water dissolved gas concentrations and fluxes can be modeled by considering an average seasonal ebullition rate in the steady state model. • Bubbles in the sediments are in local equilibrium with the dissolved pore-water gases (see Brennwald et al. 2005, section 3.1). Other than gas removal (stripping), we do not consider other physical effects of bubbles in the sediment or their release (i.e., no bubble-mediated pore-water advection; see Flury et al. 2015). We assume that bubbles leaving the sediment do not interact notably with the surrounding sediment and pore water while leaving the sediment (i.e., no additional mass transfer once they begin to exit the sediment). • Advective transport is not considered (i.e., only cohesive sediments are considered).
• Bioturbation is not explicitly considered.
• CH 4 is the only gas produced in the sediments, i.e., no N 2 production (Kipphut and Martens 1982;Martens et al. 1998;Wilkinson et al. 2015;Schmid et al. 2017). • Boundary conditions: At the SWI, the concentrations of dissolved pore-water gas are equal to the concentrations in the lake and N 2 is in equilibrium with the atmosphere. At the bottom of the sediment, a zero flux condition is applied.

Transport and reaction equations
The transport equations in our model of dissolved concentrations C i (z) of gases i in the sediment pore water are written in summarized form in Eqs. 1a and 1b. Transport processes are different for the upper nonebullitive region (Eq. 1a; sediment depth z ≤ z eb, min ) and the lower ebullitive region (Eq. 1b; z > z eb, min ; see Fig. 3): In the upper nonebullitive region, diffusive transport (first term, Fick's second law) and production (second term, W i (z)) are in equilibrium, whereas in the ebullition layer, an additional third term for the loss of gas due to bubble stripping is added. This additional stripping term is the local total ebullition rate E(z) multiplied by the individual gas mole fraction at sediment depth z. Although we only consider CH 4 and N 2 , the transport equations for any dissolved gas i can be described in the following equation.
where φ is the sediment porosity (−), D i is the effective molecular diffusion (D m ) corrected for tortuosity (D i = θ −2 D i,m ; Boudreau 1997), C i is dissolved concentration (mol m −3 ), E is total gas ebullition per bulk volume (mol m −3 d −1 ), W i the gas production per bulk volume (mol m −3 d −1 ), K H is Henry's law volatility constant (Pa m 3 mol −1 ), and P is the local critical gas  Fig. 3. Sketch of steady state sediment model based on Bazhin (2003Bazhin ( , 2010. CH 4 production W(z) leads to ebullition in the lower ebullition layer, where the total dissolved gas pressure from N 2 + CH 4 is equal to the local pressure. N 2 is continuously removed from the sediment due to stripping in the lower ebullition layer, while it is resupplied to the sediment pore water from the SWI. Gas concentrations are constant at the SWI, and fluxes are zero at the base of the sediment.
pressure (Pa). We measured porosity φ (see Supporting Information section S8 and Fig. S10) and applied a constant measured value of φ = 0.9, and therefore for tortuosity, following Boudreau (1997), a value of θ 2 = 1ln (φ 2 ) = 1.2. Temperature-dependent molecular diffusion coefficients D i,m were obtained by piecewise linear interpolation of the data presented by Jähne et al. (1987), and Henry's law volatility constants K H, i were calculated at in situ temperature according to the parameterization presented by Sander (2015). The local critical gas pressure P is calculated as P = ρgh + P atm − P H 2 O , with the water density ρ = 1000 kg m −3 , gravity g = 9.81 m s −2 , the water depth h in (m), the atmospheric pressure P atm , which is~944 hPa at Lake Soppen, and the temperature-dependent water vapor saturation pressure P H 2 O calculated with the equation of Buck (1981). The model (Eq. 1) is driven by the production of gases W i (z) (mol m −3 d −1 ). In this study, we consider no production of N 2 (W N 2 = 0), and CO 2 contributes negligible amounts to gas pressure (not modeled). We apply an exponential profile of CH 4 production rate per bulk volume (W CH4 ) as a function of sediment depth z: where a is the maximum CH 4 production (mmol CH 4 m −3 d −1 ) at the upper-most layer of sediment and b is the decay parameter (m −1 ).

Ebullition layer condition E
Within the lower ebullition layer (Fig. 3), the total dissolved gas pressure is equal to the hydrostatic plus atmospheric pressure minus water vapor saturation pressure (Eq. 3).
Bubbles are assumed to be formed only below a sediment depth z eb, min where the total dissolved gas pressure reaches local oversaturation. The two regions are described by a different set of differential equations (see Eq. 1). At the interface, concentrations and fluxes are continuous (Bazhin 2010). Before any pore-water concentrations C i (z) can be modeled, the extent of each region (i.e., z eb, min ) must be found. The depth z eb, min is defined by the combination of Eqs. 1-3. From this it follows that at z = z eb, min , the depth where bubbles start to form, the condition in Eq. 4 must be met (Bazhin 2010 the derivation of Eq. 4 is provided in Supporting Information section S9).
In Eq. 4, P res is introduced to simplify the equation and is the local atmospheric + hydrostatic pressure minus the partial pressures of gases at the SWI, which decreases with the lake water CH 4 concentration C CH4,lake : The depth z eb, min can be found numerically from Eq. 4. Note that, in general, z eb, min is dependent on the gases that are produced/consumed and on how the production is distributed (see Supporting Information section S9). In this work, Eq. 4 is valid for only CH 4 production, which decays exponentially with depth (see Eq. 2). Equation 4 is used both to determine the value of z eb, min and whether bubbles can be formed, i.e., if z eb, min is shallower than the total sediment depth (L), then bubbles are formed, otherwise bubbles are not formed.

Boundary conditions
At the bottom of the sediment z = L, a zero-flux condition is set for all gases, ∂C i ∂z z = L ð Þ= 0. The position of the lower boundary was set at L = 5 m, as we found that results did not change significantly for L > 5 m. At the SWI (z = 0), the CH 4 concentration is assumed to be equal to the concentration of the overlying lake water. The dissolved partial pressure of N 2 is assumed independent of water depth and therefore fixed at the SWI to a value of 0.78P atm (see measurements by Horn et al. [2017]).
Interestingly, the N 2 boundary conditions imply a lower limit to the amount of CH 4 in bubbles. Assuming negligible N 2 production, the dissolved gas pressure of N 2 is limited to a maximum of 0.78P atm . As bubble gas consists mainly of N 2 and CH 4 (Casper et al. 2000;Walter et al. 2008; see also our own results), the main contribution to dissolved partial gas pressure comes from these two gases, and other gases can be neglected. Considering that the required minimum total dissolved gas pressure for bubble formation increases linearly with water depth (Eq. 3; Fig. 1), then it follows that the minimum dissolved partial pressure of CH 4 required for bubble formation also increases linearly with depth. Therefore, the minimum CH 4 fraction in a sediment bubble in equilibrium with the pore water as a function of water depth h is described as:

Solution properties
Equations 1-4 are solved numerically for the concentration profile C i (z). Besides concentration profiles, the below-defined solution properties are calculated as outputs (see summary key variables Table 1).
From the combination of Eqs. 1 and 3, the total bulk gas ebullition rate with sediment depth E(z) (mmol m −3 d −1 ) is defined as Additionally, from this ebullition rate profile E(z) (Eq. 7), we can derive the sediment depth z eb,50% above which 50% of the depth-integrated ebullition occurs.
Diffusive fluxes F i,D (mmol m −2 d −1 ) of dissolved gases at the SWI are defined with Fick's first law Integrating the production rate (Eq. 2) simplifies the total production to Ð and the CH 4 ebullition pathway fraction f E of total production is Finally, the CH 4 bubble composition X CH4 is given as a ratio of the CH 4 ebullitive flux divided by total gas ebullitive flux note that F N 2 = − F N 2 ,D = F N 2 ,E , as for steady state, the bubble transport of N 2 gas out of the sediment must equal the diffusive flux back to the sediment (no N 2 production; see Eq. 1). Also note that in Eq. 12, the denominator is the total ebullition gas flux, i.e., F tot,E = F N 2 + F CH4,E . The total ebullition gas flux can also be written as F tot, E = F CH4,E =X CH4 . In the upper diffusive layer (Fig. 3), CH 4 production is balanced by loss from the sediment due to molecular diffusion. For an overview of processes, see Fig. 3.

Overview model inputs and outputs variables
A summary of the input variables that drive the model and output variables are given in Table 1 and are summarized in Fig. 3. In the general model, a and b are known inputs (Eq. 2) and the remaining variables are solved. Using an inverse modeling approach (discussed below), any two parameters from Table 1 need to be known to solve for the remaining variables (except for the combination of f E and X CH4 , which are not independent as we show in the results).
Solving for the output variables in Table 1, an example of a model solution for 5 C, 20 m water depth, and zero lake CH 4 concentration is shown in Fig. 3. Gas concentrations are constant at the SWI and fluxes are zero at the bottom end of the sediment. Ebullition E(z) of N 2 + CH 4 is fuelled by CH 4 production W(z) and occurs only in the ebullition layer where the total dissolved gas pressure of N 2 + CH 4 is equal to the local pressure. While N 2 is continuously removed from the sediment via exiting bubbles, it is diffusively resupplied from the SWI. As we assume that no significant N 2 is produced in the sediments, dissolved concentration of N 2 decrease linearly within the diffusive layer, in contrast to CH 4 , which is produced according to W(z). The model is solved with MATLAB boundary value problem solver bvp4c. Note that the model can be easily adapted to track additional gases, e.g., including a production term for N 2 or N 2 O, which may be significant in some systems (Higgins et al. 2008;Baulch et al. 2011). For code description and possible adaptions, see Supporting Information section S10 and example scripts (available via e-mail). See Supporting Information Figs. S11 and S12 for inclusion of N 2 production.

Inverse modeling approaches for CH 4 fluxes
For a given water depth and temperature, the model is driven by the generally a priori unknown two sediment production parameters a and b, which makes the model a system with two degrees of freedom. However, for any pair of observations known in Table 1, a unique pair of a and b can be found such that the two given observations are precisely matched by the model. If an assumption on either parameter a or b can be made, then only one additional observation of any of the parameters on Table 1 must be known to determine the remaining unknown parameters. Furthermore, as shown by Bazhin (2010), for any accurately measured dissolved gas sediment profile, it is possible to find the solution for a and b and all the remaining parameters listed in Table 1. We used the MATLAB global minimum search routine fmincon to find the solutions for fitting given observations (see description of code in Supporting Information section S10 and example scripts, available via e-mail).

Applied approaches to resolve CH 4 fluxes
Although there are many combinations of observations that can be investigated for resolving the CH 4 mass balance, we explore the following approaches based on our measured data. All these approaches utilize the bubble CH 4 composition (fulfilling one of the two required observations).
Approach 1: From three deployed funnel gas traps (Fig. 2) and nearby sediment gas samples, we have observations of ebullition flux at the SWI and bubble CH 4 content (described below). Given these two observations, the parameters a and b and the complete sediment CH 4 mass balance are resolved with inverse modeling.
In this approach, we assume that the gas sampled around the funnels reflects the mean ebullition flux of that funnel. To avoid disturbing the sediment just below the funnel, we sampled bubble gas in the immediate vicinity, thus the depths at which the samples were obtained slightly varied from the funnel depth. As hydrostatic pressure has a major impact on ebullition (West et al. 2016), we linearly interpolated ebullition rates between the funnels at 15 and 20 m. The depth around the 10 m funnel did not vary much, so the bubble samples were collected close to 10 m depth.
Approach 2: For the transect sediment gas samples (Fig. 2), the ebullition fluxes are not known. It is, however, possible to estimate z eb,50% from measurements of the isotopic signature δ 13 C CH4 of both the bubble gas and the sediment pore-water dissolved CH 4 . With this additional indirect observation of z eb,50% combined with the measured CH 4 bubble content, it is also possible to estimate the CH 4 mass balance using the inverse modeling.
Approach 2 is based on the idea that isotope measurements are a good method for determining the average sediment depth at which bubbles originate and thus a good estimate for either z eb,min or z eb,50% in the model. Basically, the measured bubble δ 13 C CH4 is assumed to be the same as the dissolved pore-water δ 13 C CH4 where the bubble was formed. Consequently, the corresponding isotope ratio from the sediment profile gives an approximate depth where the bubble was formed. As previously discussed, we used the averaged δ 13 C CH4 pore-water profiles from three depth ranges (Supporting Information section S7). For approach 2, we assume that these averaged profiles are representative throughout the lake (within the defined water depth ranges).
Approach 3: It is also possible to assume a value for the CH 4 production parameters a or b (Eq. 2). Here, we propose to assume either a or b as found from the application of approach 1 nearby the funnel sites and assume the selected value applies to the entire lake. Then, with additional bubble gas content data collected on the transect, we can estimate the CH 4 mass balance at the locations where the samples were collected.

Sediment area-weighted average fluxes
Finally, using these approaches, we estimate the whole-lake hypolimnetic (below 8 m) CH 4 production and sediment CH 4 fluxes. With different approaches 1-3, we obtained flux estimates of CH 4 ebullition, diffusion, and production for each sediment gas sample collected over a range of depths. From these, we can extrapolate the hypolimnetic basin-wide CH 4 flux estimates using the sediment-area-weighted averaged fluxes according to ΔA sed,i is the sediment fraction in 0.5 m intervals (Δz = 0.5 m) at depth water z i multiplied by the linearly interpolated flux F i at that depth.

Results
We present results of direct estimates of the sediment CH 4 mass balance consisting of production rates (from sediment incubations), diffusive fluxes at SWI (from pore-water profiles), and ebullition (from funnel gas trap). Furthermore, we present measured composition and sediment bubble gas and inferred sediment depth of bubble origin. Using these data, we then show the model results for the sediment CH 4 mass balance using our three approaches. We provide model results at the bubble gas sampling locations as well as integrated results for estimates of the basin-scale hypolimnetic CH 4 budget.

Result I-Observations
Net sediment CH 4 production rates Laboratory incubations of sediments were performed to estimate potential sediment CH 4 production. Figure 4 shows that net production rates decline rapidly within the first 10 cm in the sediments (see also Supporting Information Tables S1 and S2 and section S2). Following the assumption that profiles could be described by exponential decay (Eq. 2), we found parameter a ranging from 129 to 502 mmol CH 4 m −3 d −1 (mean a = 239 mmol CH 4 m −3 d −1 ) and parameter b ranging from 10.1 to 43 m −1 (mean b = 30 m −1 ; see Supporting Information Table S3). Integrated production rates corresponding to fitted a and b (Prod = a/b) range between 3.2 and 12.8 mmol CH 4 m −2 d −1 (black solid circles in Fig. 5).

Diffusive and ebullitive CH 4 fluxes at the sediment-water interface
Diffusive CH 4 fluxes at the SWI were estimated with Fick's first law (see Supporting Information section S3). The individual pore-water profiles of dissolved CH 4 and δ 13 C CH4 and estimated local diffusive fluxes at the SWI are summarized on Supporting Information Fig. S7. Resolved diffusive fluxes where bin averaged and were (mean AE 1 SD) 9.2 AE 1.7 mmol m −2 d −1 at 12 m water depth (n = 4) and 16.6 AE 11.4 mmol m −2 d −1 at 25 m water depth (n = 5; Fig. 5; Supporting Information Table S4).
Total gas ebullition fluxes (i.e., CH 4 + N 2 ) at the SWI were measured using deployed funnels gas traps at three sites (10, 15, and 21 m water depths). We corrected the data for dissolution and measured CH 4 bubble gas fraction (see Supporting Information section S4). This resulted in average estimated CH 4 ebullitive fluxes from July to December 2017 of 2.5 AE 1.0, 0.57 AE 0.12, and 1.2 AE 0.4 mmol CH 4 m −2 d −1 for the funnels at 10, 15, and 21 m, respectively (Fig. 5). For temporal dynamics, see Supporting Information section S4 and Table S5.

Composition and origin of bubble gas
We sampled sediment gas bubbles in situ with the sediment pore-water gas sampler described in the Methods section and Supporting Information section S5. Samples consisted mainly of CH 4 and N 2 , which make up on average 98.9% of the bubble gases, while CO 2 contributed only 0.77% and O 2 was present only in trace amounts (see data in Supporting Information Table S6 and section S6). Figure 6a shows that the samples at locations near the funnel sites (colored symbols) and from the transect (open symbols) all contain more CH 4 than the theoretical minimum CH 4 fraction X CH4, min defined by Eq. 6. Samples from the transect generally had higher CH 4 fractions than samples nearby funnel sites at the same water depth. Additionally, we have analyzed δ 13 C CH4 of sediment gas samples and found values in the range of −74‰ to −66.5‰ ( Fig. 6b; Supporting Information Table S6). Samples from the transect generally had higher stable isotope ratios (i.e., less negative) than samples nearby the funnel sites (Fig. 6c).
To determine depth of origin of the sampled bubble gas, we measured the δ 13 C CH4 of pore-water dissolved CH 4 (see pore-water profiles in Supporting Information Fig. S7 and section S3). Figure 6d shows the average δ 13 C CH4 for porewater profiles obtained between 8 and 13 m water depth (n = 5) and 18-26 m water depth (n = 5). We took the averages of these profiles for the depth range between 13 and 18 m. In general, the profiles of all δ 13 C CH4 were consistently lower (more negative) with sediment depth z ( Fig. 6d and Supporting Information Fig. S7). These data are used to infer the depth of origin of sampled sediment-released bubbles (Fig. 6c) by matching the bubble δ 13 C CH4 isotope values (Fig. 6b) with the corresponding δ 13 C CH4 from the pore-water profile (Fig. 6d). An example of the estimation of sediment depth of bubble origin from δ 13 C CH4 is shown in Supporting Information Fig. S9 and section S7.

Result II-Modeling sediment CH 4 fluxes Estimating CH 4 ebullition pathway
We demonstrate that, with the model, we can estimate the proportion of the CH 4 transport pathways, i.e., the relative amount of CH 4 production leaving the sediment as ebullition, with simply the knowledge of bubble gas CH 4 content (Eqs. 11 and 12; Fig. 6a). As a first step in our modeling exercise, we developed the contour (1), and (d) core 26 m (2). Assuming exponential production (Eq. 2), the best-fitting parameters a and b and the derived depth-integrated production Prod = a/b are given as value AE standard error as well as the adjusted R 2 of the exponential fit. For data, see also Supporting Information Tables S1 and S3 and Figs. S2-S5. shown on Fig. 7a. We ran the model from 0 to 26 m depth at 1 m interval varying a and b (i.e., production a/b) in Eq. 1 from a = 1-1500 mmol m −3 d −1 and b = 5-70 m −1 . The model was run using the temperature and dissolved CH 4 profiles shown on Supporting Information Fig. S1.
By running all combinations, we see that for a given depth the proportion of the total CH 4 production emitted via ebullition f E (see Eq. 11) is directly related to the bubble composition X CH4 (Eq. 12). This pathway f E has a single solution for a given bubble CH 4 fraction at a given depth. Therefore, the ebullition pathway can be estimated in every case simply with knowledge of the bubble CH 4 fraction and vice versa. As an example, Fig. 7a shows contour lines of the ebullition proportion of total production f E , which is only dependent on the water column depth and bubble CH 4 fraction X CH4 (Eq. 12) overlain with the points of measured bubble CH 4 fraction from Fig. 6a. Therefore, for a constant water depth, an increase in bubble CH 4 content corresponds to an increase in the  August 2017. The proportion of ebullition pathway relative to the total CH 4 production can be determined from bubble gas CH 4 fraction and water depth. Also shown is the minimum possible CH 4 fraction, which is dependent on water depth only and coincides with zero ebullition. (b) Ebullition pathway corresponding to measured bubble CH 4 fraction (Fig. 6a) calculated for exact individual temperatures and lake water dissolved CH 4 .
proportion of the ebullition pathway. Similarly, for a fixedbubble CH 4 content X CH4 , the ebullition pathway decreases with increasing water depth. Figure 7a also shows the theoretical minimum CH 4 fraction (Eq. 6) which coincides with the contour of zero ebullition pathway and is independent of temperature. Figure 7b shows the exact ebullition pathway proportion corresponding to measured bubble CH 4 fraction, calculated at local water temperature and dissolved CH 4 concentration.

Estimating complete sediment CH 4 balance
We have shown that it is possible to estimate the CH 4 ebullition pathway f E (Eq. 11) based on measured bubble CH 4 fraction alone X CH4 (Eq. 12). However, this approach gives no quantitative information on production or flux magnitudes. For estimates of CH 4 ebullition, diffusive flux and production, additional observations are required. To do this, we have applied the three approaches described in the Methods section. Figure 8 shows the inverse modeling results based on fitting the measured sediment bubble CH 4 content at different sampling locations in the lake and from different estimation approaches. Therefore, all the modeling results shown in Fig. 8 reproduce exactly the measured bubble gas CH 4 content previously presented in Fig. 6a. For each model calculation, the sum of CH 4 ebullition flux (Fig. 8a) and CH 4 diffusive flux at the SWI (Fig. 8b) are equal to the integrated sediment CH 4 production (Fig. 8c). The ratio of CH 4 ebullition flux (Fig. 8a) over the integrated sediment CH 4 production ( Fig. 8c) is equal to the ebullition pathway previously presented in Fig. 7b. Figure 8d shows the modeled sediment depth of bubble origin z eb,50% . Each of the approaches are detailed and compared to each other in the following section (Fig. 9).

Approach 1
For this approach, we use data from the three sites with funnel gas traps 10, 15, and 21 m, where we have measurements of ebullition fluxes at the SWI (Fig. 5) and CH 4 content in bubble gas sampled nearby the funnels (Fig. 6a). We used a linear interpolation to infer ebullitive fluxes between the 15 and 21 m funnels. Again, we assume that the ebullition flux measured at the 10 m funnel is representative for gas samples collected in the immediate vicinity of the funnel. We therefore use the inverse modeling approach with CH 4 ebullition rate and bubble gas content to infer the CH 4 production parameters a and b (Eq. 2) and thus solve for the total sediment CH 4 mass balance (i.e., production, diffusion, and ebullition) and depth of ebullition layer z eb,50% (Fig. 8). This approach is then compared to approach 2 below (Fig. 9).

Approach 2
Approach 2 is based on our idea that the sediment depth of bubble origin inferred from bubble and pore-water isotope measurements (Fig. 5c) is a good proxy for z eb,50% in the model, that is, the sediment depth above which 50% of the depth-integrated ebullition occurs. Approach 2 was applied to both sediment bubble gas samples from nearby the funnel sites and on the transect (Fig. 8). Briefly, approach 2 uses CH 4 bubble content and the isotope estimated z eb,50% to estimate parameters a and b, and thus solve the complete mass balance (Fig. 8a-d). We compare the estimates for ebullition (Fig. 9a)   Fig. 8. Modeling results (a-d) based on measured sediment bubble gas CH 4 content at different sampling locations (near funnel are solid, transects are open symbols) for approaches 1 and 2. Open red square symbols in (a-c) are measured data as in Fig. 5 and in (d), red squares were obtained from matching δ 13 C CH4 of bubbles and pore-water profiles. In (a), error bars show standard deviation of measured CH 4, in (b) standard deviation of binned diffusive fluxes, and in (c) standard error of production from exponential fit of incubations. and diffusion (Fig. 9d) from approach 1 to those obtained with approach 2 for the near-funnel data where we have measured ebullition.

Approach 3
To infer the total sediment mass balance from sediment bubble gas composition, it is possible to assume a value for either parameter a (approach 3a) or parameter b (approach 3b; Eq. 2) based on the measured sediment CH 4 production rates (Fig. 4). However, a (surface production rate) or b (depth-decay parameter) can also be estimated with approach 1 (combination of ebullition and bubble gas data). This has the advantage of fitting simpler parameters with a priori assumption of total production rate. Here, as we did not measure ebullition on the transect, we assume a value for either parameter a (approach 3a) or parameter b (approach 3b; Eq. 2) from approach 1 and apply these for the transect data.
For samples near the 15 and 21 m funnel, the estimated parameters a and b were similar to each other, whereas the 10 m funnel was somewhat different (see Supporting Information Fig. S13). For samples nearby the 15 and 21 m, we found a mean value for parameter a of 295.1 mmol m −2 d −1 and a mean value for parameter b of 27.1 m −1 . Both values are in the same order of magnitude as inferred from sediment incubations (Fig. 4), and we subsequently use those values for approach 3. In Fig. 9, we compare the ebullition and flux estimates from approaches 3a and 3b with approach 2 for the transect data.

Integrating estimated CH 4 fluxes to the lake basin
The results from the approaches 1, 2, 3a, and 3b are areaweighted averaged and integrated over the entire lake basin below 8 m water depth as basin-scale fluxes (Table 2).

Discussion
We present a modeling approach to estimate sediment CH 4 production and flux pathways based on relatively easy-tocollect field data. The basis of the approach is the fact that gas bubbles alter the pore-water dissolved gas concentrations and that the sediment bubble gas content reflects the intensity of ebullition (Chanton et al. 1989). As a reminder, our key assumptions for this approach are: (1) the system is in steady state, (2) exponential decay rates of production (although this can be changed), and (3) results represent fluxes over a seasonal scale.

Estimating CH 4 pathways from bubble CH 4 content
As shown on Fig. 6a, the theoretical minimum CH 4 content (negligible ebullition fluxes) of a bubble can be expressed as a function of depth. Based on Eq. 6, the minimum bubble CH 4 content must increase with depth and is only a function of total pressure (independent of CH 4 source but assuming no significant production of N 2 or gases other than CH 4 and bubbles in equilibrium). Approaching greater depths, the CH 4 bubble content asymptotically approaches 100% of the total bubble gas content. For example, at 100 m depth, the fluxes (d-f) in mmol m −2 d −1 from different estimation approaches 1, 2, 3a, and 3b at the different sampling locations. Symbols: Red, green, and blue are data collected near the 10, 15, and 21 m funnels, respectively, and open symbols are transect data (Fig. 2). Table 2. Overview of sediment area-weighted CH 4 flux and production estimates from different locations and estimation approaches for the hypolimnion.

Ebullition at SWI
Diffusion at SWI Integrated production Approach 1 (near funnels) 1.5 12.7 14.2 Approach 2 (near funnels) 0.8 7.9 8.7 Approach 2 (transect) 2.1 10.9 12.9 Approach 2 (all results) 1.4 9.4 10.8 Approach 3a (transect) 1.6 9.5 11.2 Approach 3b (transect) 2.0 10.2 12.1 minimum fraction of bubble CH 4 must be greater than 93%, whereas at 1000 m depth, the minimum CH 4 fraction in a bubble must be 99.2% (neglecting other potential biologically produced gases). As ebullition rates increase, however, CH 4 content is altered. For example, at a fixed depth in a lake, the higher the ebullition rate, the higher the bubble CH 4 content becomes and the greater the fraction of the total CH 4 production emitted as ebullition (Fig. 7b). Therefore, a major finding of this work is that by knowing the CH 4 bubble content and the water depth at which it was collected, the proportion between diffusive and ebullitive flux pathways can be estimated (Fig. 7b). It is also apparent from Fig. 7b that the resolution of this approach decreases with increasing depth. This means that as the depth increases, the uncertainty of our inverse modeling approaches also increases. All the collected bubble gas content data (except a single point) were in excess of the minimum CH 4 content predicted by Eq. 6 (see Figs. 6a, 7a). We can therefore infer that ebullition was ongoing at every sampling site, and at every site, we could thus quantify the fraction of production emitted as ebullition.
Using only the gas bubble composition is useful for finding local emission hotspots, lake basin-wide emission variability, and relative ebullition flux rates (Fig. 7); however, it does not give any quantitative information on the CH 4 production or fluxes. Combined with a basin-scale hypolimnetic CH 4 mass balance (e.g., Schmid et al. (2017)), the ebullition and diffusive fluxes can be easily determined. If a second parameter, together with bubble gas content, or any two parameters shown in Table 1 are determined, then the hypolimnetic lake CH 4 mass balance and flux pathways can be estimated with very good approximation (note that the combination of f E and X CH4 are not independent and need a third parameter).

Comparison of the inverse modeling approaches
Approach 1 fit the model to observations of bubble gas CH 4 content and total ebullition flux at SWI collected with in situ bubble gas traps. We found that the samples nearby the funnels 15 and 21 m with a mean value for parameter a of 295.1 mmol m −2 d −1 and a mean value for parameter b of 27.1 m −1 were comparable to the values found from sediment incubations (Fig. 4). However, the fitted parameters a and b for samples nearby the 10 m funnel fell out of this range, with parameter a above 500 mmol CH 4 m −3 d −1 and b over 40 m −1 . This can also be observed on Fig. 9a,b when comparing to the results from approach 2.
Although overall agreement with the measured data is very good (Fig. 8), we speculate that the differences observed at the 10 m site might be due to locally enhanced sediment CH 4 production rates or may be an artifact due to the inappropriate selection of physical parameters (temperature, porosity, and tortuosity) used for this method. As this site is the shallowest, it is more exposed to changing seasonal temperature and dissolved oxygen in the overlying water. Also, we cannot rule out an influence of potential N 2 production at this funnel site (for effect of N 2 production on bubble gas composition and pore-water profiles; see Supporting Information Figs. S11 and S12). Finally, a limitation of approach 1 is that it can only be applied around the areas where the ebullition data were collected.
Approach 2 is particularly interesting as all the three sediment CH 4 parameters (ebullition, diffusion, and production) are resolved entirely with independent measurements, i.e., bubble CH 4 content and bubble depth of origin. This is based on our presented hypothesis that the bubble depth of origin can be inferred with the combination of the bubble δ 13 C CH4 value with the corresponding depth of the matching δ 13 C CH4 value in the pore water ( Fig. 6b-d). This means that approach 2 can potentially be applied to quantify the CH 4 production and flux pathways without using funnel gas traps.
CH 4 production and flux estimates from approach 2 agree extremely well with data collected at all depths (Fig. 8) and also with estimates from approach 1 for depths below 10 m 9a,d). This finding also suggests that this is a robust approach to apply to the transect data. Comparing the funnel-only estimates from approach 2 with the transect-only estimates, we conclude that ebullition along the transect is up to 2-3 times higher than nearby the funnel sites, whereas diffusion about 1.2 times and production 1.3 times higher on the transect (see Table 2; Fig. 8). These results highlight both the location of relative hotspots and the need for extensive spatial coverage for accurate basin-scale estimates of CH 4 production and flux pathways.
Approaches 3a and 3b were applied using the mean values of parameters a and b (Eq. 2) obtained from modeling data (approach 1) based on samples collected around the 15 and 21 m funnel. Both approaches agreed well with approach 2 (Table 2), however, approach 3b matched better. If funnel traps are deployed, then approaches 3a and 3b only need additional data for bubble gas composition from different locations. Thus, spatial variability of fluxes and hotspots can be estimated quickly and cost effectively. If we take approach 2 as a baseline for comparison (Fig. 9), approach 3b (where we assume a basin-wide constant decay rate of production) appears to give more consistent estimates than approach 3a (where we assume basin-wide constant maximum surface production). This makes sense, as approach 3a only assumes the shape of the production curve rather than values of surface production (a) (discussed below).

Modeling basin-scale CH 4 production and flux pathways
We scale our results to provide the basin-wide modelinferred hypolimnetic CH 4 budget ( Table 2). We want to point out the differences in the estimates using only the funnels vs. the more robust estimate of including the transect data where production and fluxes were generally higher. These results certainly highlight the robustness of this approach for estimating basin-scale CH 4 production and flux pathways.

Meaning of maximum sediment production (a) and decay rate (b)
For a visualization of how CH 4 dynamics vary with depth and production, we fix the parameter b to 27.1 m −1 (i.e., approach 3b) while varying a (maximum production rate) and depth. The results of the exercise are shown as contours of bubble CH 4 content and ebullition pathway calculated as a function of depth and production ( Fig. 10; calculations based on temperature and CH 4 profiles shown in Supporting Information Fig. S1). Note, that as parameter b is fixed as a constant, any variation of integrated sediment CH 4 production corresponds to a variation in parameter a (Prod = a/b). Figure 10 clearly shows that water depth majorly influences sediment CH 4 dynamics. This also demonstrates the tendency of reduced ebullition (or no ebullition) with increasing depth (West et al. 2016).We would like to point out that our theoretical contours of ebullition pathway look strikingly similar to the empirical contours of ebullition probability presented by West et al. (2016) in their Fig. 3.
The condition for the possibility of bubble formation (zero ebullition contour on Fig. 10) can be found algebraically from Eq. 4 by setting z eb,min to infinity (Bazhin 2010). Setting z eb,min to infinity, neglecting water vapor partial pressure and solving for parameter a, we obtain the critical minimum a (maximum CH 4 production; Eq. 2) for ebullition to occur for a given decay rate b (see Eq. 2).
The minimum necessary integrated sediment CH 4 production rate for bubble formation is then consequently Prod min = a min /b. Both Prod min and a min increase linearly with water depth for constant b, temperature, and lake CH 4 concentrations. In a real lake, however, the overlaying CH 4 concentration (discussed below) and temperature influence the value of Prod min , although the linear contribution of water depth is dominant. This effect is visible in the almost linearly depth-dependent Prod min in Fig. 10. Figure 11 shows the importance of a and b for the flux pathways. Here, we provide an example of model solutions using a fixed integrated production rate, i.e., we vary a in increments of 150 mmol m −3 d −1 and change in b to maintain a constant production rate of 15 mmol m −2 d −1 . The simulation was run for 5 C, 20 m water depth, and zero lake CH 4 concentration. We see that for this constant integrated CH 4 Fig. 10. Contour plot of bubble CH 4 content (Eq. 12) and CH 4 ebullition pathway (as a proportion of total CH 4 production; Eq. 11) with a constant sediment CH 4 production decay parameter b = 27.1 m −1 (Eq. 2). As in this case, b is fixed as a constant and variation in the integrated sediment CH 4 production corresponds to a variation in parameter a (Prod = a/b). Note also that the minimum CH 4 production necessary to sustain ebullition for a given parameter b increases approximately linearly with water depth. Fig. 11. Example sediment pore-water model calculations for 5 C at 20 m water depth. Three cases (a-c) with varying parameters a and b but same total CH 4 production. production, the more rapidly the CH 4 production rate decays with sediment depth (i.e., high parameter b), the lower the maximum pore-water CH 4 concentration becomes, despite a higher a. Consequently, the CH 4 fraction in bubble gas and the CH 4 ebullition pathway are reduced as a and b increase.
Simplistically stated, for constant sedimentation rates and organic matter content, the parameters a and b in an exponential decaying sediment CH 4 production profile (Eq. 2) can be explained by the remineralization of the more labile organic matter over time due to methanogenesis following first-order kinetics. The parameters a and b would then reflect both the kinetics of methanogenesis and the sedimentation rate including the organic matter content and quality. We propose the relationship for sedimentation with constant organic matter input (quality) a / k 1 × C l, ini / k 1 × R PP /w, b / k 1 /w, and Prod = a/b / R PP , where k 1 is a first-order reaction constant, C l, ini reflects the initial concentration (at the SWI) of labile carbon source to be converted to CH 4 , R PP is the primary production, w is the sedimentation rate, and Prod is the integrated sediment CH 4 production. Therefore, we come to a paradoxical conclusion that for a constant sedimentation rate and organic carbon input, the more refractory the carbon is (low k 1 within certain limits), the higher the CH 4 ebullition will likely be. Similarly, but not as counterintuitive, increasing the sedimentation rate (high w) of a constant labile carbon source would also lead to increased ebullition rates. This is in line with observations reported by Sobek et al. (2012) and Maeck et al. (2013). Of course other factors also affect ebullition, including depth, oxygen exposure, and temperature (Sobek et al. 2012).

Effect of overlying CH 4 water concentration on flux pathways
It is worth noting here that the analyses performed in this article used constant boundary conditions (temperature and overlying CH 4 concentration). CH 4 production is a function of temperature, and will therefore change when the overlying waters warm over the stratification season. This is especially important in shallow littoral systems where temperatures can increase substantially. Also important for the flux pathwaysfraction emitted as either CH 4 diffusive or ebullitive flux-is the overlying CH 4 concentrations. As shown in Fig. 12 and assuming constant CH 4 production using the example from Fig. 11b, as the overlying water dissolved CH 4 concentration increases over time, the fraction emitted as diffusive fluxes will decrease and the ebullitive fraction will increase. As diffusive fluxes are a function of the concentration gradient at the SWI, increasing the overlying CH 4 concentration will decrease the concentration driving force, and thus the diffusive flux. Assuming a constant CH 4 production, this translates into a higher ebullitive flux. Thus, a remediation process, such as hypolimnetic oxygenation (McGinnis et al. 2004) could reduce ebullition by simply maintaining low hypolimnetic CH 4 concentrations.

Hotspots and intra-lake variability of CH 4 fluxes
Hotspots in terms of the relative value of CH 4 ebullition pathway can be detected based solely on the composition of sediment bubble gas as shown in Fig. 7. To quantify the CH 4 ebullition flux, however, one of the presented approaches must be applied to estimate the whole sediment CH 4 balance. This also helps to quantify the spatial variability of ebullition and diffusive CH 4 flux within a lake. To obtain a complete picture of basin-scale variability of CH 4 ebullition and diffusive fluxes, we propose to use a combination of several of the presented approaches when feasible. Any of the approaches combined with data collected at several point locations (e.g., bubble gas funnels) increases the robustness of basinscale flux estimates.

Key assumptions and applicability to other waterbodies
The proposed model is applicable to other freshwaters where the same key assumptions apply: (1) the model conceptualization (Figs. 1, 3), (2) cohesive sediments (i.e., no advection), (3) negligible bioturbation, (4) exponential CH 4 production, (5) no N 2 production and N 2 concentration at the SWI is at atmospheric equilibrium. We discuss each of these assumptions below.
1. The model used in this study has been tested for validity both in the laboratory (Kusmin et al. 2006) and with a modeling study based on literature data (Bazhin 2010). Therefore, we assume that the basic conceptualization (two distinct layers) and impact on dissolved gas is applicable to many freshwater sediments. In its current form, it is a steady state model; therefore, the model represents averaged bubbling events (season scale) and lacks the resolution to resolve very short time scales. Fig. 12. Example calculation to show the effect of increasing CH 4 concentration in overlying water column on CH 4 ebullitive flux and fraction of production emitted as ebullition (ebullition transport pathway). Example calculated for same parameters a and b, temperature, and depth used for simulation in Fig. 11b.
2. The sediments are cohesive and we assume no advection. Advective transport arising from externally impressed flow (e.g., groundwater and near-bottom currents) is generally not present for muddy aquatic sediments because of their low permeabilities (Boudreau 1997). Advective transport arising from sedimentation and compaction can be neglected for the considered time and length scales (see e.g., Strassmann et al. 2005). The model is potentially applicable to permeable sediments with some modifications (McGinnis et al. 2014), however, more work is needed. 3. Bioturbation is not addressed explicitly, however, can be included simply by increasing the effective diffusion coefficient in the bioturbated zones. The inclusion of other types of bioturbation may need more modeling effort (Boudreau 1997). 4. We have applied an exponentially decaying CH 4 production assuming anoxic conditions. Numerically, it is straight forward to change the production profile W(z) (Eq. 2) for any algebraic or numerical function (see Supporting Information section S10). For systems with oxic overlying water, we expect oxygen penetration in the top millimeters of the sediments (Bryant et al. 2010). We do not model CH 4 oxidation explicitly; however, the model formulation (Eq. 1) is still valid under these conditions. A possible formulation was presented by Wilkinson et al. (2015), who considered combined SWI effects with a sigmoid function. Such a general W(z) could even be negative (net consumption) in the top parts. Finally, the laboratory results of Kusmin et al. (2006) and the model of Bazhin (2010) suggest that the model is applicable under oxic conditions. 5. We assume that N 2 production is negligible for the interpretation of our bubble gas composition samples below 13 m from June to December 2017. For the samples above 13 m, we cannot rule out the possibility of some influence due to denitrification. In lakes with high N 2 production in the sediments (for example, due to denitrification, anaerobic ammonium oxidation, or denitrifying anaerobic methane oxidation), both the approaches and the model used to estimate CH 4 emissions applied in this study would require additional parameterizations. The model itself can be modified to include N 2 production by including an N 2 production term or profile. As an example, we show model simulations using an exponentially decaying sediment N 2 production (Supporting Information Figs. S11 and S12). This adaption of our model would be applicable for the study of sediments where both denitrification and methanogenesis take place (Higgins et al. 2008).

Implications for CH 4 flux measurement methods
With approach 1, it is possible to close the sediment CH 4 mass balance at sites in the lake where funnel gas traps are deployed without any additional pore-water measurements; hence, approach 1 is already an improvement of the funnel gas trap method alone. Referring to Table 1 with possible combinations for estimation approaches, we can see that we can combine a measurement of total ebullition rate with, e.g., the measurement of diffusive flux CH 4 at SWI or sediment CH 4 production. This applies not only to the funnel gas trap method but also to hydroacoustic methods. Both methods rely for the calculation of CH 4 ebullition flux on the fraction of CH 4 in bubble gas, which is a priori not known but can be resolved for with a combination of any two parameters (Table 1).

Summary
We have confirmed the validity of measurement-modeling approaches with which it is possible to resolve the complete sediment CH 4 flux mass balance and basin-scale production and flux pathways. Measuring bubble gas composition alone already allows for resolving the relative proportion of the flux pathways at a given depth (i.e., diffusive versus ebullitive CH 4 fluxes) and identification of ebullition hotspots. Combining sediment bubble gas composition with additional measurements of ebullition flux at SWI (approach 1), depth of sediment bubble origin based on δ 13 C CH4 signature (approach 2), or sediment CH 4 production parameters a and b calibrated with approach 1 (approach 3a/3b) allows for solving the complete mass balance of CH 4 and flux pathways.
To summarize, approach 1 allows to close the sediment CH 4 mass balance at sites in the lake where funnel gas traps are deployed, without any additional pore-water measurements. Conversely, approach 2 is a novel way to estimate the entire sediment CH 4 mass balance, entirely independent from funnel gas trap measurements. Approach 2 is applicable to the whole lake but does require the measurement of δ 13 C CH4 signature in pore waters and bubble gas. Finally, if funnel gas traps are installed, approach 3 allows to quickly assess the CH 4 fluxes of the rest of the lake in cost-effective manner. All of these measurement approaches will allow better quantification of CH 4 ebullition and diffusion at SWI and integrated sediment CH 4 production.
For hypolimnetic CH 4 budgets based on evolution of CH 4 concentrations and knowledge of basin-scale diffusivity, Schmid et al. (2017) demonstrated that while the production rate of CH 4 can be estimated, the fractional contribution of ebullition could not be elucidated based on these data alone. Performing such a budget with the addition of sediment bubble gas composition, thus adding estimates of the flux pathways, would allow closing the CH 4 budget with higher certainty. In the case of Schmid et al. (2017), this could be relatively easily achieved by collecting data on sediment bubble gas composition at the various locations over time, following the recommendation for sampling frequency and locations of Wik et al. (2016). Such additional data would also improve the mass balance by accounting for spatial heterogeneity, and allow for accurate estimates of CH 4 bypassing the water column and emitted to the atmosphere via ebullition.
Current ebullitive CH 4 emission estimates from inland freshwaters are highly uncertain given their stochastic nature. Applying our approaches, combined with hypolimnetic CH 4 mass balances, will certainly allow the identification of CH 4 hotspots, better resolve in-lake CH 4 budgets, and improve overall CH 4 emission estimates from inland freshwaters and their contribution to global CH 4 budgets.