Role of gas ebullition in the methane budget of a deep subtropical lake: What can we learn from process‐based modeling?

We analyzed the processes affecting the methane (CH4) budget in Lake Kinneret, a deep subtropical lake, using a suite of three models: (1) a bubble model to determine the fate of CH4 bubbles released from the sediment; (2) the one‐dimensional physical lake model Simstrat to calculate the mixing dynamics; and (3) a biogeochemical model implemented in Aquasim to quantify the CH4 sources and sinks. The key pathways modeled include diffusive and bubble release of CH4 from the sediment, aerobic CH4 oxidation, and atmospheric gas exchange. The temporal and spatial dynamics of dissolved CH4 concentrations observed in the lake during 3 years could be well represented by the combined models. Remarkably, the relative contributions of ebullition and diffusive transport to the accumulation of CH4 in the hypolimnion during the stratified period could not be accurately constrained based only on the observed evolution of CH4 concentrations in the water column. Importantly, however, our analysis showed that most (∼99%) of the CH4 supplied to the water column by bubble dissolution and diffusive transport from the sediment is aerobically oxidized, whereas a substantial fraction (∼60%) of the sediment‐released bubble CH4 is directly transported to the atmosphere. Ebullition is thus responsible for the bulk of the emissions from Lake Kinneret to the atmosphere. Therefore, as in all freshwaters, ebullition quantification is crucial for accurately assessing CH4 emissions to the atmosphere. This task remains challenging due to high spatio‐temporal variability, but combining in situ measurements with a process‐based modeling can help to better constrain flux estimates.


Introduction
Methane (CH 4 ) concentrations in the atmosphere have increased by a factor of 2.5 since pre-industrial times (Dlugokencky et al. 2011). Between 1999 and 2006 concentrations remained nearly constant, and it was thought that a new steady state had been reached (Heimann 2011). However, concentrations have been rising again since 2006 (Nisbet et al. 2014). The reasons for these dynamics are not well understood and have been controversially debated (Aydin et al. 2011;Bousquet et al. 2011;Kai et al. 2011;Levin et al. 2012;Simpson et al. 2012). Significant uncertainties remain in the quantification of the sources and sinks of atmospheric CH 4 (Kirschke et al. 2013;Ghosh et al. 2015) and further efforts are required to better estimate both natural and anthropogenic sources of CH 4 as well as to investigate feedbacks between global change and CH 4 emissions. CH 4 emissions from lakes and reservoirs contribute significantly to the global CH 4 sources to the atmosphere (Barros et al. 2011;Bastviken et al. 2011). CH 4 is produced under anoxic conditions in freshwater sediments during the degradation of organic matter either from acetate fermentation or CO 2 reduction (Conrad 2005). It can then be emitted to the atmosphere via different pathways (Bastviken et al. 2004): formation and subsequent release of CH 4 bubbles from the sediments (ebullition); diffusive transport from the sediment to the hypolimnetic water, where CH 4 can accumulate under anoxic conditions and then be transported to the surface by turbulent mixing or to downstream rivers by deepwater withdrawal from reservoirs (Gu erin et al. 2006); and plantmediated emissions (Carmichael et al. 2014).
The major process limiting CH 4 emissions from lakes is microbially mediated aerobic or anaerobic CH 4 oxidation in the sediment ( a Norði et al. 2013) or the water column (Hakemian and Rosenzweig 2007;Caldwell et al. 2008). Methane oxidation often effectively constrains diffusive emissions from deep stratified lakes, where only a small fraction of the CH 4 released to the hypolimnia of these lakes reaches the atmosphere (Schubert et al. 2010).
Ebullition fluxes, which have been highlighted as a major pathway for CH 4 emissions (DelSontro et al. 2010;Grinham et al. 2011), are spatially heterogeneous and intermittent and therefore difficult to quantify. The spatial variability of gaseous CH 4 fluxes closely corresponds to the highly variable gas content in the upper sediment layer that, in turn, depends on organic matter content and granulometric characteristics of sediments (Katsnelson et al. in press). CH 4 bubble formation in the sediment occurs when production is faster than CH 4 oxidation and diffusion (Katsman et al. 2013). These bubbles can eventually overcome the sediment tension and be released to the water column (Boudreau 2012). A fraction of the emitted gaseous CH 4 is subsequently dissolved in the water column, the amount of which depends on the bubble release depth, the bubble size and dissolved gas concentrations in the surrounding water, while the rest is directly emitted to the atmosphere (McGinnis et al. 2006). The release of bubbles from the sediment strongly depends on the hydrostatic pressure (Joyce and Jewell 2003). Ebullition fluxes can therefore vary by orders of magnitude depending on whether the lake level is rising or falling (Ostrovsky 2003;Scandella et al. 2011;Ostrovsky et al. 2013).
The aim of the present study is to quantify the main sources and sinks of CH 4 and its delivery to the atmosphere from a medium-size subtropical lake. We chose Lake Kinneret as a case study, because it has characteristics of both a lake and a reservoir and its intensive monitoring provides the necessary information on physical driving forces and lake chemistry for modeling. Using state-of-the-art modeling, we assess the role of CH 4 ebullition from the sediment in the whole-lake CH 4 budget and its emission to the atmosphere. Since the water level of Lake Kinneret is strongly influenced by water inflow and withdrawals and has varied by up to 5.5 m over the last decades (Rimmer et al. 2011a), we also investigate the effects of changing water levels on ebullition. We further analyze the possibility to discriminate between the diffusive and ebullitive fluxes based on the evolution of CH 4 concentrations in the water column. Finally, we conclude with a summary of the implications and uncertainties of the model approach and suggestions to better constrain the ebullition contribution to the water column CH 4 budget, and to atmospheric emissions.

Study site
Lake Kinneret is a warm, deep monomictic lake located in the northern part of the Afro-Syrian Rift Valley (Fig. 1). The 2730 km 2 Lake Kinneret watershed includes four hydrological units with large differences in typical precipitation. Lake Kinneret is one of the most studied deep subtropical freshwater lakes in the world . A summary of the present scientific knowledge on the lake has been published in a recent book .
The water and salt balances of Lake Kinneret have been studied in great detail (Rimmer and Gal 2003). The major inflow to the lake is the Jordan River, which drains $1700 km 2 of the watershed. At its maximum water level (-209 m a.s.l.) the lake's volume is 4200 3 10 6 m 3 , its surface area is 167 km 2 , and its maximum depth is 44 m. The main outflow from the lake is via Israel's National Water Carrier, withdrawing 100-400 3 10 6 m 3 of water annually. Imbalances between natural inflows and outflow of water have resulted in amplified water level fluctuations in the last few decades which are expected to further increase due to growing water demand and climate change (Rimmer et al. 2011b).
The lake is thermally stratified from April to December, with surface water temperatures exceeding 308C in recent years (Imberger and Marti 2014). The epilimnion thickness varies seasonally from $10 m in spring to 20-25 m in autumn, and complete mixing usually occurs by the end of December (Rimmer et al. 2011a). Anoxic conditions in the hypolimnion start to develop in the bottom boundary layer and in the lower metalimnion in May, subsequently spreading over the whole hypolimnion and prevailing throughout the stratified period (Nishri et al. 2000(Nishri et al. , 2011. At the same time, phosphate and reduced substances such as ammonium, sulfide and CH 4 accumulate in the hypolimnion (Nishri et al. 2000;Eckert and Conrad 2007). Lake Kinneret is mesoeutrophic with a mean annual net primary production of $600 g C m 22 (Yacobi 2006;Ostrovsky et al. 2013).
Due to the high energy content in the internal wave field, horizontal currents are fast, resulting in time scales for horizontal mixing in the hypolimnion of only a few days (Marti and Imberger 2008). Except for depth-dependent boundary processes, such as sediment resuspension and redistribution (Ostrovsky and Yacobi 2010), biogeochemical processes in the hypolimnion can thus be assumed to be horizontally homogenous (Ostrovsky and Yacobi 2009). One-dimensional models can therefore be used to simulate these processes in the lake.
Previous modeling studies have addressed various physical and biogeochemical processes in Lake Kinneret (Pan et al. 2002;Gal et al. 2003Gal et al. , 2009Bruce et al. 2006;Rimmer et al. 2006a;Vernieres et al. 2006;Shilo et al. 2007;G omez-Giraldo et al. 2008;Li et al. 2013). However, none of these studies assessed the sources and sinks of CH 4 in the lake.

Data
The models developed in the following section utilize data collected by various research groups, and made available in the Lake Kinneret Database . For driving the models, meteorological data, lake surface temperatures, the water and salt balance of the lake, as well as light absorption data in the epilimnion are required. Observations of temperature, CH 4 and oxygen (O 2 ) concentrations, and CH 4 ebullition are used for parameterization and calibration. In the following, we shortly summarize the sources of these data. More detailed information, especially concerning the meteorological data and the water and salt balance, are given in Section S1 of the Supporting Information.

Meteorological data
Air temperature, wind speed, and wind direction from the Tabgha Station ) were used for driving the physical lake model (Fig. S1 in the Supporting Information). The station is located $8 m above the lake surface, 700 m offshore (Fig. 1). Since the surface boundary condition for temperature is forced to observed temperatures rather than calculated from heat fluxes, no other meteorological data, such as humidity, cloud cover, or longwave radiation, were required for the simulations.

Lake water chemistry and temperature
The following data used in this study were all measured at the deepest point of the lake at station A ( Fig. 1) as part of the routine monitoring of the Kinneret Limnological Laboratory (KLL) and retrieved from the Lake Kinneret Database . These data include for the period between January 1998 and February 2001: 157 vertical profiles of O 2 , measured using the Winkler method (APHA 1971); 73 observations of light attenuation in the surface layer of the lake (Yacobi 2006); 143 vertical profiles of conductivity and temperature (CTD profiles), measured using a STD-12 Plus (Applied Microsystems). Lake surface temperature was interpolated in time from the average values in the top 1 m of the temperature profiles. In addition to the routine monitoring, 64 vertical profiles of dissolved CH 4 concentrations were measured between June 1998 and February 2001 with triplicate water samples followed by headspace analysis using a Shimadzu GC with an FID detector (for the dataset and details on the CH 4 measurement method see Eckert and Conrad (2007)).

Water and salt balance
A detailed monthly water and chloride balance of Lake Kinneret was presented by Rimmer and Gal (2003) and Rimmer and Givati (2014). The respective in-and outflows and the inflow salinities calculated from chloride concentrations were considered in the physical model, but neglected in the methane model, as their contribution to the overall CH 4 budget of the lake is small (see "Discussion" section). The total water balance was negative (-110 3 10 6 m 3 yr 21 ) during the simulated period, resulting in a lake level drop of 2.5 m (Fig. 2). More details on the water and salt balance are given in Section S1 of the Supporting Information.

Bubbles
To study the depth-dependence of bubble density near the bottom, eight acoustic surveys were carried out in (11 September 2008, 24 September 2008, 05 August 2009, 09 September 2009, 14 October 2009, 25 November 2009, 04 August 2010 September 2010) using a Biosonics dual-beam Scientific Echosounder DE5000 (6.5 8 half-power narrow beam width) operating at 120 kHz, a pulse width of 0.2 ms, and a ping rate of 5 pings s 21 . The acoustic data were collected along two transects connecting the central station A with the northern tip of the lake and with the MS location, respectively (Fig. 1). Acoustic data collected in the anoxic hypolimnion were echo-integrated to estimate the volumetric density of bubbles in the absence of fish (Ostrovsky 2009a). The echo integration procedure was also used for the calculation of bubble density at shallower locations where fish and bubble targets coexist. The relative contributions of bubbles and fish echoes were quantified based on a separation of bubble tracks from fish tracks using data collected with a slow boat speed (Ostrovsky 2009b).
For each transect, the volume backscattering coefficient, s v , was calculated in bins of 3-m width (0.5-3.5 m above the bottom) and 500 pings (100 s sampling duration). The volumetric density of targets per cubic meter (N) in each bin was estimated by scaling s v with the mean backscattering crosssection r bs , i.e., N 5 s v /r bs (Ostrovsky 2003). Each observation was referred to the mean bottom depth of the bin. The overall mean N was calculated for 1-m bottom depth intervals for all surveys. The depth-specific mean N were computed for all surveys and were based on 4-19 samples (bins). The ratio of standard deviation to the mean was usually close to 1, reflecting significant variability of bubble density over the sampling period. The total number of samples used for these calculations was 338. The bubble flux (F, no. of bubbles m 22 s 21 ) emanating from the sediment surface was calculated as F 5 Nv, where v is the mean rise velocity of bubbles of 0.22 m s 21 (Ostrovsky 2003).
Bubble size distributions were determined from acoustic surveys during five sampling campaigns (09 December 2010, 07 July 2011, 11 August 2011, 14 September 2011 November 2011) using a 120 kHz split-beam Simrad EK60 scientific echosounder (78 beam width transducer). The data were collected in the near-bottom 3-m stratum positioned in the anoxic hypolimnion where no fish are present. The sampling scheme was the same as described above. Target strengths (TS 5 10 log(r bs )) of 2634 bubble echoes collected in 14 transects ( Fig. 1) were converted into bubble volume using an empirically determined TS-volume relationship (DelSontro et al. 2015).

Model introduction
The CH 4 balance of Lake Kinneret was simulated using a combination of three different models. Vertical profiles of CH 4 concentrations were simulated with a basic onedimensional reaction-diffusion model using the software AQUASIM 2.1g (Reichert 1994). This model requires vertical diffusivities as a function of time, which were calculated externally with the mixing model Simstrat (Goudsmit et al. 2002). The dissolution in the water column of CH 4 bubbles emanating from the lake sediments was also calculated externally using a discrete bubble model (McGinnis and Little 2002;McGinnis et al. 2006). These three models are described in the following sections, and are called methane model, physical model and bubble model, respectively. An overview of the links between the three models (Supporting Information Fig. S2) and the availability of the model code (Supporting Information section S3) are presented in the Supporting Information.
The overall modeling approach used here is comparable to the bLake4Me model developed by  for simulating CH 4 emissions from shallow pan-arctic lakes . In comparison, our model includes a more sophisticated physical transport model and considers the effects of lake bathymetry, with sediment surface areas and thus CH 4 sources simulated at all depths within the lake. In the gas transport module of bLake4Me, the lake area is constant with depth and gas sources occur only in the lowest water column grid cell. In return, bLake4Me includes a module that explicitly simulates the thermal dynamics, as well as CH 4 production and bubble formation in the sediment. Including such a module in a model with depthdependent bathymetry is much more complex and computationally expensive, as a separate instance of the sediment module has to be run with different boundary conditions at each depth within the lake, as it was done in the model LAKE 2.0 recently developed by Stepanenko et al. (2016). We therefore decided to use a simpler approach for simulating the CH 4 sources from the sediment. The difference in the model setups between bLake4Me and our model can be explained by the different systems investigated in these studies. Whereas for a deep stratified lake such as Lake Kinneret, a correct representation of the lake-internal vertical structures is indispensable, the simulation of CH 4 emissions from mostly shallow and seasonally frozen pan-arctic lakes requires a stronger focus on the processes within the sediments.

Physical model
A modified version of the Simstrat model was used to calculate coefficients of vertical diffusivity based on the budget equations for turbulent kinetic energy (k), representing the source of energy for turbulent mixing, the energy dissipation rate (e) and the vertical stratification. Besides the classical k sources, shear and buoyancy, Goudsmit et al. (2002) introduced an additional term to account for boundary mixing by internal seiche motions, which improves the ability of the model to reproduce turbulent diffusivity in the hypolimnia of lakes. This model has been used to successfully reproduce mixing in qualitatively different deep lakes such as Lake Baikal , Lake Ohrid (Matzinger et al. 2007), or Lake Geneva (Perroud et al. 2009), and also performed reasonably well without calibration for lakes in different climate zones in the lake model intercomparison project LakeMIP (Stepanenko et al. 2010(Stepanenko et al. , 2013(Stepanenko et al. , 2014Thiery et al. 2014). The modifications compared to the version of Simstrat described by Goudsmit et al. (2002) are presented in section S4 of the Supporting Information.
Simulations with the physical model were run from 15 February 1998 to 24 February 2001. The model calculations were performed for each year separately, with the model being initialized with the vertical temperature profiles as observed on 15 February 1998, 22 February 1999, and 14 February 2000. February was chosen as the time of initialization because the lake is usually almost completely homogenized at this time of the year. Salinity variations in the vertical profiles were negligible on these dates, and the average salinity measured was used as initial condition for the whole water column. The lake area was defined as a function of depth based on data by Ben-Avraham et al. (1990).
All inflows (described in section S1 of the Supporting Information) were introduced into the system between 24 m and 28 m, while the outflows were removed between 27 m and 25 m. Thus, both inflows and outflows are below the lake surface level at all times during the simulations. Since the epilimnion is well-mixed during most of the year down to these levels, the results are insensitive to the exact input and output depths. The inflows also introduce heat and salt which were added to the model based on existing heat and salt balance calculations (Rimmer et al. , 2011b. The heat input is treated in the model as a temperature input and was determined by multiplying the total inflows from rivers and rain with the temperature of the river Jordan, and the inflows from groundwater sources with an average temperature of 268C. Salt was introduced with the inflows and modified by evaporation and mineralization of organic material as described in section S1 of the Supporting Information. The parameters of Simstrat were set to the default values as described by Goudsmit et al. (2002) with the following exceptions: three parameters were used as fit parameters during the calibration procedure (see "Results" section); the stability functions as described by Galperin et al. (1988) were switched on, which improved the simulation of the mixing depth during the stratified period.

Bubble model
The bubble model was used to calculate the source term for dissolved CH 4 in the methane model that results from the dissolution of bubbles during their ascent in the water column, and the ebullition of CH 4 to the atmosphere at the lake surface. It includes the release of bubbles from the lake sediments, which varies with both depth and time, and the exchange of CH 4 and other gases (i.e., O 2 and N 2 ) between the bubbles and the water phase.
The source term in the methane model (see Eq. 5) due to CH 4 dissolved from bubbles, r CH4,b (mmol m 23 s 21 ) at a certain depth z is calculated by integration over depth of the CH 4 dissolved from bubbles released at all depths below z, using the following equation: (1) Here, F CH4,b (mmol m 22 s 21 ) is the bubble flux of CH 4 from the sediment averaged over the entire lake area, which is a calibration parameter in the methane model, f t,CH4 (-) is a time-dependent function of the lake-level variation (step 4 below), f z,CH4 (-) is a depth-dependent function of bubble release (step 2 below), A surf is the lake surface area and A the lake area as a function of depth, f d,rel (mm 21 ) is the probability density function (pdf) of the contributions of each bubble diameter to the total CH 4 release from the sediment (step 1 below), f d,diss (z 0 , z rel , d) is the fraction of CH 4 dissolved at a given depth z 0 as a function of the release depth z rel , and the initial bubble diameter d (mm) (step 3 below). The integrals were approximated by summing up over 0.2-m depth intervals and 1-mm steps for bubble diameters between 3 mm (d min ) and 20 mm (d max ).
The source term r CH4,b is calculated in four steps: 1. The bubble size distribution ( Fig. 3) was determined from the average distributions observed during five measurement campaigns as described in the data section. It is assumed to be constant with both time and depth. It has been previously shown that using just one bubble size with the Sauter diameter should give similar results as simulating the entire bubble size distribution (McGinnis and Little 2002;DelSontro et al. 2015). However, the bubble size distribution would still have to be known in order to calculate the Sauter diameter. 2. The bubble flux per area from the sediment is derived from measurements of bubble size, rise velocity, and the volumetric density of bubbles near the bottom at different depths (see "Data" section). It varies as a function of depth due to depth-specific hydrostatic pressure change, variability of sediment properties, organic matter accumulation rate, and CH 4 production rate. It is multiplied by the sediment area available within each depth interval and normalized to yield the pdf of the bubble source as a function of depth (Fig. 4). 3. The dissolution in the water column of CH 4 bubbles emanating from the lake sediments was determined using a discrete bubble model (McGinnis and Little 2002;McGinnis et al. 2006). The model tracks individual bubbles rising in water and accounts for gas dissolution and stripping of nitrogen, O 2, and CH 4 based on background and bubble concentrations. The change in bubble size with depth is thus continually updated and reflects the change in mass of gas within the bubble, changing water temperature, and decreasing hydrostatic pressure. The new bubble size is then used to calculate the bubble rise velocity and mass transfer coefficient. The dissolution rate of single pure CH 4 bubbles is calculated for bubble diameters ranging from 3 mm to 20 mm with 1-mm intervals released at depths between 5 m and 45 m depth with 5-m intervals. The vertical distribution of the CH 4 dissolution is not sensitive to the initial fraction of CH 4 in the bubbles. For example, assuming an initial composition of 90% CH 4 (Ostrovsky et al. 2008) and 10% N 2 instead of 100% CH 4 , the total gas flux (and thus the number of bubbles with the same bubble size distribution) would simply have to be increased by 11% to achieve the same flux of CH 4 , resulting in a very similar vertical distribution of the CH 4 release. Background concentrations of 100 lmol L 21 CH 4 below 20 m depth and 0 lmol L 21 above were used, but the results of the model are not sensitive to these assumptions within the range of observed concentrations. The concentrations of the other dissolved gases in the water column have no significant effect on the dissolution rate of the bubbles within the concentration ranges observed in Lake Kinneret, therefore constant concentrations were used in the model calculations. Dissolved nitrogen was assumed to be at equilibrium with the atmosphere. Background temperatures and O 2 were obtained from the Lake Kinneret Database. 4. The bubble source term varies with time due to changes in lake level, as a lower pressure at the sediment surface promotes the formation and release of CH 4 bubbles. Hydroacoustic measurements (Ostrovsky et al. 2008) showed that massive ebullition of CH 4 took place in 2001 and 2002, when the lake water level was the lowest ever recorded. A rapid rise in lake level in winter-spring 2002-2003 and a corresponding increase in hydrostatic pressure stopped gas emission from the sedimentary gas reservoirs in 2003. Increasing free gas content in the sediment with decreasing water levels were confirmed by evaluating acoustic properties of the sediment (Ostrovsky and Tegowski 2010). The function f t,CH4, which describes the timedependence of the bubble flux from the sediment was calculated as a function of time based on the water level history. The regression equation previously derived by Ostrovsky et al. (2013) was used to calculate f t,CH4 as a function of the difference (DLL) between the current lake level and the average lake level in the months July to November during the previous year with the exponential factor k t,CH4 5 20.8247 (6 0.068, r 2 5 0.97). The regression is based on measurements in 6 yr, spanning DLL between 21.5 m and 13.5 m and bubble densities ranging over two orders of magnitudes. f t,CH4 is normalized with the constant c n to average to unity during the years 1998 to 2001. The response of the model to k t,CH4 is tested in a sensitivity analysis where this parameter is increased or reduced by 50%, i.e., far outside the uncertainty of the regression.

Methane model
Methane concentrations and fluxes were calculated with a one-dimensional diffusion-reaction model using the lake module of the software AQUASIM 2.1g (Reichert 1994). This software is designed to simulate physical mixing and biogeochemical processes in natural and artificial aquatic systems. A vertical grid of 0.25 m resolution was used. The time step was adjusted by the software depending on stability criteria for the solution of the differential equation system. The model numerically solves the equations: where [O 2 ] and [CH 4 ] (mmol m 23 ) are the concentrations of dissolved oxygen and methane; A (m 2 ) is the lake crosssectional area as a function of depth; K Z (m 2 s 21 ) the vertical diffusivity calculated as a function of time and depth with the physical model, and smoothed using a running average with a window of 5 d in order to increase computation speed and numerical stability (Fig. S3 in the Supporting Information); z (m) is the vertical dimension positive upwards; r CH4 and r O2 are the source/sink terms for CH 4 and O 2 . The source/sink term r O2 includes both O 2 production by primary production and O 2 consumption in the water column. O 2 production was simulated with a constant rate r O2,prod determined by model calibration in the surface layer down to a depth of 10 m. O 2 consumption was simulated as a constant rate r O2,cons also determined by model calibration and limited to a maximum first order consumption with a rate of 10 25 s 21 (corresponding to an exponential decay with a time scale of 1 d) in order to avoid negative concentrations.
The source/sink term for CH 4 is given by: The sources include r CH4,b which accounts for the CH 4 dissolved from rising bubbles (Eq. 1), and r CH4,d , the diffusive flux of dissolved CH 4 from the sediment surface. The latter is given by: for all depths where the O 2 concentration is below a threshold value of 0.1 mg L 21 , and set to zero otherwise. F CH4,d is the areal diffusive flux of CH 4 out of the sediment and is a model fit parameter. The third term in Eq. 5 represents aerobic CH 4 oxidation, which is implemented with Michaelis-Menten kinetics. v max,CH4 is the maximum CH 4 oxidation rate, and K M,CH4 and K M,O2 are the half-saturation constants for CH 4 and O 2 . Most literature values reported for K M,CH4 are in a relatively narrow range of 2-10 lmol L 21 (Hanson and Hanson 1996;Knief and Dunfield 2005), whereas v max,CH4 can vary over several orders of magnitudes, with reported values ranging at least from 0.3 lmol L 21 d 21 for water samples from an arctic lake (Lofton et al. 2014) to > 1000 lmol L 21 d 21 in some samples from the sedimentwater interface of Lake Washington (Lidstrom and Somers 1984). Since the two parameters have similar effects on CH 4 oxidation rates at the concentrations observed in the epilimnion of Lake Kinneret, only one of them can be used as a calibration parameter. We therefore fixed K M,CH4 at 5 lmol L 21 , and used v max,CH4 as a calibration parameter. K M,O2 has  (dotted), of the release of methane in bubbles from the sediments (thick dashed, calculated from the product of the sediment area at each depth and the depth-dependent bubble flux) and of the dissolution of methane from bubbles rising in the water column. The integral over depth of the pdf for methane dissolution is only $0.4, as the remaining 0.6 is emitted directly to the atmosphere. In addition, the pdf of methane dissolution resulting from a depth-independent bubble release from the sediment is shown (thin dashed), as it was used in the sensitivity analysis.
only little effect on the simulated CH 4 concentrations (Fig. S5 in Supporting Information) and was set to a literature value of 0.7 mg L 21 (Lidstrom and Somers 1984). Anaerobic oxidation of methane (AOM) potentially contributes to the CH 4 budget of Lake Kinneret. AOM has been shown to be coupled to sulfate reduction, denitrification, and the reduction of iron or manganese (Cui et al. 2015). The processes and kinetics of AOM in the water column of lakes are largely unknown, but previous studies suggested that it is generally slow even if potential electron acceptors are readily available (Durisch-Kaiser et al. 2005;Morana et al. 2015). As a consequence, it seems unlikely that AOM removes a notable fraction of the CH 4 that has been released from the sediments of Lake Kinneret. Because of the unavailability of data to constrain the process and its likely small contribution to the overall CH 4 budget (see also "Discussion" section), we chose to neglect AOM in our model. Lightdependent aerobic CH 4 oxidation is another process that can consume CH 4 in anoxic waters . However, in Lake Kinneret this process is apparently irrelevant due to extremely low light availability below the oxycline.
Both CH 4 and O 2 are exchanged with the atmosphere at the lake surface. The flux of the species C from the atmosphere to the lake, F C (positive numbers mean flux into the lake) was quantified using the relationship: The equilibrium concentration C eq of CH 4 was set to 0.003 mmol m 23 . Its temperature dependence could be neglected since the lake surface water is oversaturated by several orders of magnitude at any time. For O 2 , C eq was calculated using Henry's law from its partial pressure in the atmosphere and its solubility as a function of temperature. The gas exchange velocity v gas (m s 21 ) was calculated as an empirical function of wind speed (Cole and Caraco 1998): v gas 5f vgas;T 5:75 Á 10 26 15:97 Á 10 27 Á u wind 1:7 where u wind (m s 21 ) is the wind speed measured at the KLL station, and f vgas,T is a factor that accounts for the effect of temperature on the Schmidt number and correspondingly on gas exchange.  (Read et al. 2012) and differ between gases due to gas-specific bubblemediated transfer (Goddijn-Murphy et al. 2015). However, the sensitivity analysis in the Supporting Information as well as the results of a test simulation with doubled gas transfer velocity showed that the choice of the gas exchange velocity does not affect any of the main conclusions of this study. One disadvantage of the lake module in AQUASIM 2.1g is that it does not allow varying the lake surface level. However, because the lake surface layer is usually well-mixed, this does not have significant effects on the results of the simulations, and a constant lake surface level at 2212 m a.s.l. was used for the simulations. The initial concentration of CH 4 (after seasonal mixing) was set to 0.1 lmol L 21 , which is a typical value observed at this time of the year. O 2 was initialized with the vertical profile observed on 15 February 1998. The temperature at the lake surface, used only for calculating gas exchange, was forced to the observed surface temperature.
Three different parameter estimations were iterated multiple times until the parameter values converged, using the parameter estimation algorithm of Aquasim in the following sequence: (1) the rates of O 2 consumption and production were fitted to minimize the root-mean-square deviation (RMSD) between observed and simulated O 2 concentrations at all depths; (2) the diffusive and bubble fluxes of CH 4 were fitted to minimize the RMSD between observed and simulated CH 4 concentrations below 20 m depth; (3) the maximum aerobic CH 4 oxidation rate was fitted to minimize the RMSD between observed and simulated CH 4 concentrations in the top 15 m of the water column. The choice of the parameters used for model calibration was based on a sensitivity analysis (section S6 in the Supporting Information).

Physical model
Three model parameters were used to fit the simulated temperatures to observations by minimizing the RMSD between observed and simulated temperatures: C 10 determines the efficiency of energy transfer from the wind to the lake. It has a default value of 0.001 in Simstrat; the factor a determines the transfer of energy from the wind to the internal seiche (default 0.01, however, usually lower); the parameter q N (default 0.75) determines how the seiche energy is dissipated vertically as a function of the density stratification. Only observations at depths > 20 m were used for the fitting, because our main fit target was to quantify the vertical exchange between the surface layer (where temperature is forced to match observations at the surface) and the deep water. Because relatively small deviations of the metalimnion depth result in large temperature deviations, fitting to observations in the entire depth range would mainly aim at reproducing the correct metalimnion depth. The parameter estimation was done using the patternsearch algorithm of MATLAB 7.10, which is a direct search algorithm that looks for a minimum of the objective function based on an adaptive mesh aligned with the coordinate directions in the parameter space.
The variation of the parameters for the three simulated years is small, indicating that the model calibrated for one year has good predictive capabilities for the other years (Table 1). Interestingly, the values of a are much lower than for other lakes where this model has previously been applied (Goudsmit et al. 2002;Finger et al. 2007;Perroud et al. 2009), suggesting a comparably lower efficiency of wind energy transfer to seiche motions and boundary mixing. Conversely, the calibrated values of C 10 were similar to those expected for an average wind speed of 2.7 m s 21 (0.0014; W€ uest and Lorke 2003).
The agreement between observed and simulated temperatures is good (Fig. 5). The mean absolute difference for the entire dataset is 0.418C, the RMSD is 0.828C, the bias 20.088C. For the deep water below 20 m depth, which was used for calibrating the model, the mean absolute error is 0.318C, the RMSD is 0.528C, and the bias 20.118C. In general, the model tends to slightly underestimate mixing and thus the warming of the hypolimnion in spring. One possible reason for this underestimation may be the disregard of the indirect warming of the hypolimnion by heat exchange with the sediment during internal seiching, which has recently been reported to contribute to spring hypolimnion warming in Lake Kinneret (Nishri et al. 2015). The increase in salinity observed during the study period resulting from the decreasing water levels (Rimmer and Nishri 2014) was replicated by the model (Fig. S4 in Supporting Information), indicating that the salt balance of the lake was correctly implemented.

Bubble model
The bubble model resulted in the vertical distribution of CH 4 release to the water column from bubble dissolution displayed in Fig. 4. It was estimated that only $7% of the CH 4 released from bubbles in the entire lake is dissolved in the water column below 20 m depth (Fig. 4), 33% is dissolved above 20 m depth, and the remaining 60% are released to the atmosphere. The functional dependence of the bubble flux on lake level variations developed by Ostrovsky et al. (2013) results in a large seasonal and inter-annual variability (Fig. 6). Estimated fluxes are typically reduced to

Schmid et al. CH 4 budget of a deep subtropical lake
half of the long-term average during rising water levels in spring, and rise up to double the long-term average at the end of sinking water levels in December. Due to the sinking water level trend, bubble fluxes were above the long-term average during the three study years. The model predicts almost no bubble fluxes for the year 2003 when the lake level was steeply rising.

Methane model
As will be discussed further below, the simulations with the methane model showed that the relative contributions of ebullition and diffusive transport to the accumulation of CH 4 in the hypolimnion during the stratified period could not be accurately constrained. For this reason, we present first the results of the "optimized version" of the methane model where both the diffusive and the bubble flux were determined by parameter estimation, and subsequently discuss the results of simulations covering a large range of F CH4,d , the model parameter determining the bubble release from the sediments.
The methane model was able to reproduce the observed O 2 concentrations in the lake (Fig. 7) using two calibration parameters. The zero-order O 2 consumption rate of O 2 in the water column was estimated at 0.081 mg L 21 d 21 . The second parameter, the O 2 production in the epilimnion, was estimated at 0.30 mg L 21 d 21 . Simulated O 2 concentrations in the epilimnion closely followed observations, except that spring peak concentrations were underestimated. The model reproduced the O 2 concentration dynamics in the deep water well, which is an additional indication for the accuracy of the simulated hypolimnetic mixing processes by the physical model. The difference most relevant for the simulation of CH 4 concentrations was that the model underestimated the decrease in O 2 concentrations at 30 m depth and below in summer 1999, which then caused a delay in the build-up of CH 4 concentrations.
Three parameters were estimated to minimize the RMSD between observed and simulated CH 4 concentrations in different depth ranges: The first parameter, v max,CH4 , was estimated to be $16 lmol L 21 d 21 by minimizing the RMSD to the observed CH 4 concentrations in the top 15 m of the water column. The model captures the average epilimnetic CH 4 concentrations, but misses some of the observed peaks during the stratified period (Fig. 8). The two other parameters were used to minimize the RMSD between observed and simulated CH 4 concentrations below 20 m depth: the average flux of CH 4 bubbles (corresponding to a relative CH 4 flux of 1 in Fig. 6) from the sediments was estimated at 34 mmol m 22 d 21 or 12 mol m 22 yr 21 ; the diffusive flux of CH 4 from the sediments was estimated at 9.5 mmol m 22 d 21 or 3.4 mol m 22 yr 21 , however, only when O 2 concentrations are < 0.1 mg L 21 at a given depth. The sediment area was exposed to O 2 concentrations < 0.1 mg L 21 during 55% of the time below 20 m depth, and during 32% of the time averaged over the entire lake. Therefore, the mean annual diffusive flux is 1.9 mol m 22 yr 21 below 20 m depth and 1.1 mol m 22 yr 21 averaged over the entire lake.
A sensitivity analysis showed that the capability of the model to independently estimate the two CH 4 source parameters is limited based on the observed data, even though the effects of the two source terms on the depth distribution of the CH 4 concentrations in the lake are rather different (Fig.  9). Increased bubble fluxes more strongly affect concentrations in the upper hypolimnion, whereas increased diffusive fluxes lead to higher concentrations in the lower hypolimnion. Higher bubble fluxes also result in higher CH 4 concentrations in the epilimnion, but this can be compensated in the model by increasing v max,CH4 .

Schmid et al. CH 4 budget of a deep subtropical lake
In order to illustrate the behavior of the model, the results of additional model runs are presented in Fig. 10. In these simulations, the bubble flux was modified between 0 mmol m 22 d 21 and 65 mmol m 22 d 21 . The model was re-calibrated for each bubble flux using the other four calibration parameters. Reasonable agreement with the observed O 2 and CH 4 concentrations in the water column was achieved for a broad range of bubble fluxes by modifying the diffusive flux source to achieve a similar total flux to the hypolimnion. Agreement with observations in the epilimnion is achieved by modifying the O 2 source and sink terms as well as the maximum CH 4 oxidation rate. When the bubble flux exceeds 65 mmol m 22 d 21 , modification of these three parameters is not sufficient anymore to reproduce the observed CH 4 concentrations in the lower part of the epilimnion.
The optimal values of the bubble and diffusive fluxes furthermore depend significantly on some of the model assumptions, as shown by a few additional simulations. The estimated bubble flux increases with increasing temporal variation of the bubble source term, as shown by two simulations where k t,CH4 in Eq. 2 was increased and decreased by 50%, respectively. The estimated bubble flux also increases with the gas exchange velocity, as shown by a simulation where v gas was increased by a factor of two, and it decreases if the bubble flux is assumed to have no depth-dependence, i.e., if the dotted line in Fig. 4 is replaced by a constant value, resulting in a different depth distribution of CH 4 dissolution from bubbles (Fig. 4). The latter is the only scenario that significantly deviates from the approximately linear dependence between the two flux terms, where a lower total flux is required without the depth dependence because a larger fraction of the CH 4 released with bubbles is dissolved in the water column.

Evaluation of the modeling approach
The comparison between observed and simulated concentrations of O 2 and CH 4 (Figs. 7,8) shows that their dynamics in Lake Kinneret can be largely explained using a restricted set of processes. For O 2 , these include a constant production in the epilimnion (down to 10 m depth) by photosynthesis, and a constant zero-order consumption in the entire water column. For CH 4 , the modeled processes are the release of CH 4 bubbles from the sediment and their dissolution in the water column, and diffusive fluxes of dissolved CH 4 from the sediment when the overlying water is anoxic. Additionally, both species are affected by aerobic CH 4 oxidation (which also consumes O 2 ) and by gas exchange across the water surface.
It was not necessary to include AOM as an additional process, which supports the conclusion of Eckert and Conrad (2007) that AOM likely plays a minor role in the water column of Lake Kinneret. This agrees with the findings from other stratified lakes, where AOM in the water column was detected but mostly assessed to play a minor role in the overall CH 4 budget (Lopes et al. 2011;Pasche et al. 2011). An exceptional case is Lake Tanganyika where the occurrence of a permanent water layer with sufficiently high sulfate and CH 4 concentrations seems to provide ideal conditions for AOM, and a significant fraction of CH 4 is therefore removed by this process despite a slow reaction rate (Durisch-Kaiser et al. 2011). AOM is generally much more relevant, however, for constraining diffusive CH 4 fluxes from the sediment to the water column a Norði et al. 2013;Deutzmann et al. 2014). In the water column of Lake Kinneret, AOM likely occurs in the layers below the oxycline, where both CH 4 and potential electron acceptors are available. The role of AOM in the CH 4 Fig. 10. Estimated model parameters (left) and root mean square differences (RMSD) between simulated and observed methane and oxygen concentrations (right) as a function of the bubble flux from the sediment. The model parameters were estimated separately for each depicted run. The colored dots indicate the following simulations: 1 (brown): optimal solution; 2 (green) and 3 (pink): simulations where the temporal variability of bubble fluxes is decreased and increased, respectively; 4 (orange) simulation with doubled gas exchange velocity; 5 (violet) simulation without depth dependent gas bubble flux. For the simulations 2-5, the optimization procedure was repeated after changing the model parameters.
dynamics could be to reduce CH 4 concentrations just below the oxycline and thus the flux of CH 4 by turbulent diffusion from the hypolimnion to the epilimnion. The quantification of this process would require profiles of CH 4 concentrations and its carbon isotopic composition with high vertical resolution in the relevant depth range, which are not available for Lake Kinneret.
The qualitatively most important deficiency of our relatively basic biogeochemical model is that it cannot capture the O 2 and CH 4 peaks that occur in spring in the epilimnion. These O 2 peaks are likely caused by peak primary production in the spring months . Similarly, local methanogenesis in the epilimnion linked to primary production (Grossart et al. 2011;Tang et al. 2016), a process that is not considered in our model, is one hypothesis to explain the observed shallow CH 4 peaks, with a further possible contribution by lateral transport from littoral zones (Hofmann et al. 2010). In contrast to the fast lateral homogenization in the hypolimnion, the time scale for lateral mixing in the epilimnion of Lake Kinneret is on the order of 2 months (Imberger and Marti 2014), allowing the formation of horizontal patches of phytoplankton (Ng et al. 2011). Similarly, high local release of CH 4 near the boundaries may cause significant lateral variation in epilimnetic CH 4 concentrations, which cannot be considered in a one-dimensional model.

Sources of methane in Lake Kinneret
A surprising outcome of this study is that it was not possible to accurately quantify the relative importance of the diffusive and ebullition fluxes in Lake Kinneret only based on observed vertical profiles and temporal dynamics of CH 4 concentrations. As shown by the sensitivity analysis above, a reduction in one of these two fluxes can be compensated in the model by increasing the other flux, without creating notable changes in the vertical patterns of CH 4 concentrations in the water column. This has important consequences for the greenhouse gas budget of the lake, as a large fraction of the ebullition flux is emitted to the atmosphere while most of the diffusive flux is oxidized in the water column. To account for this uncertainty, we first discuss the best solution resulting from the model optimization procedure and then discuss the CH 4 budget of the lake for a wide range of the two CH 4 sources.
The best overall agreement with the observed dissolved CH 4 and O 2 concentrations was achieved with an annual ebullition flux between 15 February 1998 and 15 February 2001 of 12.3 mol m 22 yr 21 and a diffusive flux from the sediments of 1.9 mol m 22 yr 21 below 20 m depth (1.1 mol m 22 yr 21 averaged over the entire lake). As discussed in the following, both values are above expectations from previous assessments. Nevertheless, since the total CH 4 flux to the hypolimnion is well constrained by the observed concentrations, it is not possible that both fluxes were overestimated by the model. This is corroborated by the results of a previous study with a simpler model considering only the diffusive flux as a source term and neglecting ebullition, where diffusive fluxes at the sediment-water interface of 30 mmol m 22 d 21 (11 mol m 22 yr 21 ) and 16 mmol m 22 d 21 (6 mol m 22 yr 21 ) had to be postulated for the years 1999 and 2000, respectively (Rimmer et al. 2006b).
The diffusive flux in Lake Kinneret is about 2-4 times the values of diffusive fluxes of CH 4 measured from the sediments of some eutrophic Swiss lakes (M€ uller et al. 2012), but higher fluxes of $3.6 mol m 2 yr 21 were recently observed in small eutrophic Soppensee (Switzerland) (D. Vachon et al. pers. comm.). Diffusive fluxes above those observed in Switzerland could be expected, as mineralization rates increase with temperature (Gudasz et al. 2010), and hypolimnion temperatures in Lake Kinneret typically exceed those in Swiss lakes by $108C. The estimated diffusive fluxes are also about a factor of 4 higher than the average flux estimated from the gradient of CH 4 concentrations in sediment cores from the deepest location of Lake Kinneret Sivan et al. 2011). Additional measurements are required to evaluate whether the flux estimates at this location are representative for the entire lake, as CH 4 profiles at the sediment-water interface and thus the estimated diffusive fluxes can show a large spatial and temporal variability Steinsberger et al. 2017). Furthermore, it has been shown that the presence of gas voids can increase the diffusive flux of CH 4 in sediments as the voids provide shortcuts for the transport of weakly soluble gases . For example, the presence of 5% gas voids can increase the flux of CH 4 by one order of magnitude. Recent acoustic observations indicate important spatial variability in void fractions in the surface sediments of Lake Kinneret (Katsnelson et al. 2017). Further investigations are required to assess their effects on the average gas fluxes from the sediment.
The seasonally averaged (June-December) bubble fluxes calculated from the observed bubble densities as described by Ostrovsky et al. (2008) (Fig. 7), but they are still a factor two to three below the flux estimated by the optimized model. Furthermore, the contribution of ebullition to the total fluxes from the sediments during anoxic conditions in the hypolimnion ($80%) exceeds a previous estimate of 50-75% based on CH 4 budgets within sediment cores .
The CH 4 released from the sediment is produced from decaying organic matter, supplied by organic carbon exported from the epilimnion. The average annual carbon export from the epilimnion to the hypolimnion has been estimated to $20 mol m 22 yr 21 (Eckert and Conrad 2007;Ostrovsky and Yacobi 2010), and the observed export in the years 1998-2000 was 13-14 mol m 22 yr 21 (I. Ostrovsky, unpubl. data). Between 10% and 40% of this carbon is permanently buried in the sediment , and in the absence of additional sources of electron donors (such as H 2 ), a maximum of half of the remainder can be converted to CH 4 by methanogenesis. The estimated carbon export is thus barely sufficient to supply enough carbon for the two CH 4 sources in the model.
Several causes might contribute to these discrepancies. Due to the sinking water levels, the CH 4 fluxes were likely above the long-term average during the simulated period. Furthermore, recent observations indicate that the lake level effect, which was estimated from long term annual observations as included in the model (Fig. 6), might not fully replicate the amplitude of seasonal variations, with gas fluxes during the winter months being significantly lower than predicted by the model (I. Ostrovsky, unpubl. data). Also, the carbon export measured with sediment traps at the central station A above the benthic boundary layer (Ostrovsky and Yacobi 2010) does not account for focusing of organic particulate material via the turbulent benthic boundary layer (Ostrovsky and Yacobi 1999), which increases actual carbon deposition in the deeper parts of the lake (Ostrovsky et al. 2014). In fact, the sediment mass accumulation rate at station A has been shown to be about twice that observed at shallower locations except the area near the Jordan River inflow .
For the CH 4 budget of the hypolimnion, it makes a large difference whether the same amount of CH 4 is released from the sediment as bubbles or by diffusive transport. Only about 7% of the bubble flux is dissolved in the water column below 20 m depth (Fig. 4), whereas most of the diffusive flux accumulates in the hypolimnion during the stratified period. Consequently, if we reduce the diffusive flux in the model by, e.g., 1 mmol m 22 yr 21 , the bubble flux needs to be increased by about 10 mmol m 22 yr 21 to produce similar concentrations in the hypolimnion (Fig. 10). Increasing the contribution of the diffusive flux in the model thus implies a much smaller total CH 4 production in the sediment. The sensitivity analysis shows that the RMSD between observed and simulated concentrations of CH 4 and O 2 does not increase dramatically within a wide range of the relative contributions of ebullition and diffusion to the total flux. Since the model cannot well constrain the contributions of the two sources, and the available data on the carbon budget of the lake suggests that the total CH 4 production in the sediment should be smaller than predicted by the model, we conclude that the bubble flux estimated by the optimized version of our model is likely an upper limit.
The methane budget of the lake Figure 11 shows the calculated CH 4 budget for the hypolimnion (below 20 m depth), for the epilimnion (above 20 m depth) and for the entire lake for a wide range of bubble fluxes from the sediment. The depth of 20 m was selected to separate the epilimnion and the hypolimnion because the water contains O 2 down to this depth during most of the year. The real epilimnion thickness varies seasonally from $10 m in spring to 20-25 m in autumn (Rimmer et al. 2011a). The CH 4 budget for the optimized model solution (the brown dot in Fig. 10) is given in Table  2, including separate budgets for two seasons with different properties due to O 2 availability in the hypolimnion. From January to May, the hypolimnion is mostly oxic, whereas it is anoxic most of the time at most depths from June to December.
In-and outflows were not included in the methane model, but their contributions to the CH 4 budget ( Table 2) are negligible. Water is withdrawn from the epilimnion of the lake, where CH 4 concentrations are on average < 1 lmol L 21 . Given a total water output of $3 3 10 8 m 3 yr 21 , the CH 4 output via the outflows is $10 5 mol yr 21 . No observations of CH 4 concentrations in the inflows are available. But riverine CH 4 concentrations usually do not exceed a few lmol L 21 and tend to increase with the fraction of wetlands in the catchment McGinnis et al. 2016), which is low for the Jordan River. The CH 4 input with the inflows to Lake Kinneret is therefore unlikely to exceed 10 6 mol yr 21 .
For most of the bubble flux range, both ebullition and diffusive fluxes contribute to a similar extent to the accumulation of CH 4 in the hypolimnion during the stratified period, where the diffusive flux is much more influential during the anoxic period. During the holomixis in winter, O 2 is supplied to the sediment surface, and thus aerobic oxidation of CH 4 in the topmost sediment effectively constrains the diffusive flux to the water column from January to May . During this time, the CH 4 flux is limited to ebullition, which in Lake Kinneret is generally low during winter and spring due to the seasonal rise of water levels and slower organic matter decomposition at comparably low temperatures. The model might actually even underestimate the seasonality of the CH 4 fluxes, as it does not consider the temperature dependence of methanogenesis. Methanogenesis in ecosystems has been shown to increase with temperature with an activation energy of 0.86-1.07 eV (Yvon-Durocher et al. 2014). Assuming that this is also valid for Lake Kinneret, the resulting seasonal variations of methanogenesis in the lower hypolimnion (seasonal temperature variations of 2-38C; Fig. 5) would be $30-50%, and correspondingly higher in the shallower parts of the lake. Since the summer fluxes are well constrained by the observations, adding a seasonal temperature dependence would reduce the model-estimated fluxes during winter. The depth dependence of bubble release due to temperature decreasing with depth in summer is, however, already included in the pdf of bubble release in Fig. 4. It should also be noted that the locally adapted population of methanogens in the lake may react differently to temperature, and that a recent study has shown no increase in CH 4 concentrations in sediments of two temperate lakes with increasing temperature (Fuchs et al. 2016).
More than half of the CH 4 released to the hypolimnion is oxidized below 20 m depth during times when sufficient O 2 is supplied by seasonal mixing, whereas a smaller fraction is transported upward to above 20 m depth, mostly during the seasonal mixing in autumn and winter. Release from bubbles is the largest CH 4 source for the epilimnion, except if the bubble flux is < 6 mmol m 22 d 21 , when the input by diffusion from the hypolimnion becomes more relevant. Diffusion from the sediment is negligible in the epilimnion, as the CH 4 produced in the deeper sediments is oxidized within the surface sediment layers. Based on the model results, the dissolved CH 4 in the epilimnion is almost entirely oxidized, with only a small amount of $10-15 mmol yr 21 released to the atmosphere by gas exchange. This corresponds to 1-1.5% of the CH 4 dissolved in the entire lake during the year. This value is rather independent of the model assumptions, as long as the model reproduces the observed concentrations in the epilimnion, which hardly exceeded 1 lmol L 21 , even during the early mixing period in December (Fig. 8). However, the diffusive flux to the atmosphere increases if we assume a higher gas exchange velocity, and it has also been suggested that this flux could be enhanced by up to a factor of two or three by the formation of microbubbles in the lake surface layer (Prairie and Del Giorgio 2013;McGinnis et al. 2015), a process that is not included in our model. Nevertheless, the fraction of CH 4 dissolved in the water column that is finally released to the atmosphere likely does not exceed the low percentage range. This confirms previous findings of efficient CH 4 removal by aerobic oxidation (Gu erin and Abril 2007;Schubert et al. 2010;Zigah et al. 2015), even though it has been suggested that an important fraction of stored CH 4 can be released during the seasonal overturn in shallower lakes (Encinas Fern andez et al. 2014).
Despite the small gas exchange from the epilimnion to the atmosphere, Lake Kinneret is a large source of CH 4 due to ebullition. This is true even if the total bubble flux should be significantly lower than estimated by the optimized model. Less than 10% of the CH 4 released from bubbles are dissolved in the hypolimnion (below 20 m depth), 30% are dissolved higher above, and 60% are directly emitted to the atmosphere. This is comparable to previous findings for the Iron Gate reservoir, where most of the bubble-released CH 4 was estimated to be directly emitted (McGinnis et al. 2006). These values depend on the bubble size distribution (Fig. 3) and on the vertical source distribution (Fig. 4). But independent of the model assumptions, it is clear that an important fraction of the CH 4 released as bubbles from the sediment reaches the atmosphere. In the optimized model, about 99% of the gas release to the atmosphere stems from ebullition (Table 2). Even if this is probably an upper limit, as discussed above, it highlights the importance of ebullition as a source of CH 4 from stratified lakes to the atmosphere. Ebullition is apparently the major CH 4 source even from lakes where the rate of bubble release from the sediments is an order of magnitude smaller than estimated here for Lake Kinneret. This is in line with the conclusions of Bastviken et al. (2011), who calculated global CH 4 ebullition to be about a factor 5 higher than diffusive fluxes based on published observations, but estimated that this factor was likely higher in reality due to biased sampling. Averaged over the lake surface area, the estimated ebullition flux from Lake Kinneret corresponds to CH 4 emissions of $20 mmol m 22 d 21 . This is above the range of previous observations of ebullition from lakes of 0.02-12 mmol m 22 d 21 ), but well within the range of observations from reservoirs that reach up to $45 mmol m 22 d 21 (Deemer et al. 2016).
In total, the model results indicate that Lake Kinneret releases approximately the same amount of CH 4 by ebullition to the atmosphere as it removes by aerobic oxidation. An accurate quantification of ebullition fluxes is therefore an important prerequisite for estimating greenhouse gas fluxes not only from shallow lakes (Sepulveda-Jauregui et al. 2015) but also from deeper lakes. This is a difficult task, however, as these fluxes can vary by orders of magnitude both in space and time Bussmann et al. 2013).

Influence of lake level variations and climate change
The lake level variations in Lake Kinneret cause more than order-of-magnitude inter-annual variation of CH 4 Diffusive flux to the atmosphere -15 (-3/-12) -15 (-3/-12) ebullition (Fig. 6). Based on recent observations, the seasonal variations within one year are probably high as well. This has several relevant consequences. First of all, for a serious quantification of the CH 4 fluxes from the lake, ebullition at the water surface must be measured during different periods, but with a focus on a high temporal and spatial resolution during sinking water levels when emissions are expected to be highest. While it is clear that water level variations strongly affect the temporal distribution of CH 4 ebullition as the reduced pressure facilitates bubble formation and release (Ostrovsky 2003;Maeck et al. 2014), it still remains unknown to what extent they modify the long-term average ebullition rate by changing the preferential pathway of CH 4 out of the sediment from diffusion to ebullition. Intuitively one would expect a larger average ebullition under frequently changing pressure conditions that should lead to abrupt ejections of gaseous CH 4 from sediments. In contrast, constant hydrostatic pressure may lead to a dominance of the diffusive process, where fluxes from the sediment to the water column can be constrained by CH 4 oxidation (e.g., AOM). Still, this needs to be confirmed in laboratory and field measurements. If this is the case, lake level variations would increase overall CH 4 emissions to the atmosphere. This effect might even be amplified in the future, as higher water level fluctuations may be expected in arid and semiarid regions due to predicted possible decreases in precipitation and increasing frequency of extreme climatic events (Suppan et al. 2008;Samuels et al. 2009;Rimmer et al. 2011b) as well as stronger anthropogenic pressure on the existing water resources . Furthermore, CH 4 production and ebullition in organic-rich sediments increase with water temperature (Eugster et al. 2011;Wik et al. 2014). Thus, global warming, as well as anthropogenic factors causing increase in near-bottom temperature (e.g., discharge of hypolimnetic water by hydropower stations, and hypolimnetic aeration) may intensify CH 4 ebullition from lakes and reservoirs.

Implications, uncertainties, and future directions
The models used in the present study for the assessment of the CH 4 budget of Lake Kinneret contain specific uncertainties. For the physical model, the most significant challenges are to accurately quantify the vertical turbulent diffusivity in the hypolimnion and the depth of the metalimnion. For deep stratified lakes this can usually be achieved with sufficient accuracy if the required meteorological forcing and observed temperature and salinity profiles for model calibration are available. The transport of bubbles through the water column and the exchange of gases between bubbles and the water can be simulated with good accuracy. The major uncertainties therefore lie in the quantification of the sinks and sources of CH 4 and their spatial and temporal variation.
In stratified lakes with an anoxic hypolimnion such as Lake Kinneret, the net accumulation of CH 4 in the hypolimnion during the stratified period can be determined from profiles of dissolved CH 4 concentrations. In the absence of significant anaerobic CH 4 oxidation in the water column, this accumulation is governed by the diffusive flux from the sediment and by dissolution from rising CH 4 bubbles. It was one goal of the present model approach to assign the total flux to the two sources, based on the different vertical distribution in the water column of CH 4 dissolved from bubbles or diffused from the sediment. Surprisingly, this was not possible for the case of Lake Kinneret, as similar agreement with observations could be achieved for a wide range of relative contributions of the two sources. This implies that at least one of the two sources needs to be quantified independently, such that the other source can be estimated with the model.
Separating these two sources is crucial for estimating the CH 4 emissions from a lake, because a much larger fraction of the CH 4 released by ebullition from the sediment reaches the atmosphere. For this reason, the difficulty of the model to discriminate between the sources directly translates to a high uncertainty in estimating the CH 4 emissions to the atmosphere. Additional constraints by observed carbon fluxes within the lake indicate that the contribution of ebullition, and thus the emission to the atmosphere, was likely overestimated by the optimized model. Conversely, this apparent overestimation may also indicate that the available budget of organic carbon for Lake Kinneret could underrate the amount of organic carbon burial due to uncertainties associated with the spatial heterogeneity of sedimentation processes (e.g., Yacobi 1999, 2010). In any case, it is a robust result of the simulated CH 4 budget that ebullition is the dominant pathway of CH 4 release to the atmosphere from Lake Kinneret, contributing more than 90% to the total emissions. The high contribution of ebullition is favored by the lake level variations, and was probably above average during the investigated years due to sinking water levels.
A detected discrepancy between CH 4 accumulation in the hypolimnion of Lake Kinneret and the direct observations of diffusive and ebullitive fluxes implies that available measurements likely misrepresent the actual fluxes. Further efforts are required for an accurate experimental quantification of these fluxes. The diffusive fluxes from the sediment may be better constrained by additional measurements of porewater concentrations in sediment cores. However, these fluxes have been shown to vary spatially, e.g., due to variable sediment deposition and focusing (Carignan and Lean 1991;Steinsberger et al. 2017), as well as seasonally and interannually due to variable oxygen availability at the sediment surface and inter-annual variations of carbon deposition . A large number of porewater profile measurements at different depths and seasons are thus required to reliably quantify average diffusive CH 4 fluxes. This might be facilitated by novel measurement techniques such as porewater sampling by capillary electrophoresis (Torres et al. 2013), but especially in bubbling sediments, the determination of diffusive fluxes remains a difficult task Tyroller et al. 2016).
The spatial and temporal variability of ebullition is even higher than that of the diffusive flux. Wik et al. (2016) estimated for three small sub-arctic lakes that measurements with inverted funnels at $10 locations and on $30-40 sampling days would be required to quantify annual ebullition with an accuracy of $20%. For the much larger Lake Kinneret, even more locations would be needed for the same accuracy, and the procedure would have to be repeated over several years to assess the additional effects of lake level variations. The further development of advanced methods for monitoring and quantification of gas ebullition, using acoustic techniques (e.g., scientific echosounder, multi-beam sonar), and automated bubble traps (Wilkinson et al. 2015) and flux chambers (Duc et al. 2013), might help to more accurately determine the spatio-temporal variability of this flux. Eddy-covariance (Eugster et al. 2011;Podgrajsek et al. 2016) and hyperspectral methods (Gålfalk et al. 2017) can be used to quantify the flux from a lake to the atmosphere over larger spatial scales. Finally, the total carbon fluxes from the sediment could be better constrained by long-term measurements of the carbon deposition with sediment traps in combination with measurements of spatial variation of carbon burial in the sediments.
A basic deficiency of the one-dimensional model approach used in this study is that it cannot explicitly account for the spatial heterogeneity of the processes. This is not critical for simulating the accumulation of dissolved CH 4 in the hypolimnion, as the horizontal transport processes are fast compared to the residence time of the CH 4 (Marti and Imberger 2008). But spatial heterogeneity is relevant in the epilimnion, where both CH 4 oxidation and exchange with the atmosphere occur on time scales of days. A one-dimensional model, which is usually calibrated with measurements at the center of a lake, cannot account for local processes such as ebullition hot spots near inflows, or increased CH 4 production in the littoral zone and the resulting transport of CH 4 toward the pelagic zones and to the atmosphere. As a consequence, the assessment of the contribution of local sources in peripheral areas to the net CH 4 emissions with a one-dimensional model remains a challenging task. If observations from a sufficiently large number of lakes show clear functional dependencies between lake properties and the contributions of peripheral areas to the overall fluxes, these spatial processes may be parameterized as a function of lake properties in the future.
The presented results show that a combination of a physical mixing model, a bubble model, and a relatively simple geochemical model for the CH 4 sources and sinks is a useful tool to investigate the roles and relevance of the processes defining the CH 4 budget of stratified lakes. The exchange of information between the three models can be tedious, however, and for future investigations, it would be preferable to merge these three models into a single piece of software. Once such a combined model has been set up and calibrated for a specific lake, this approach can be expanded to project the response of CH 4 fluxes and emissions to future climatic changes, extreme events or management actions.