LAKE 2 . 0 : a model for temperature , methane , carbon dioxide and oxygen dynamics in lakes

A one-dimensional (1-D) model for an enclosed basin (lake) is presented, which reproduces temperature, horizontal velocities, oxygen, carbon dioxide and methane in the basin. All prognostic variables are treated in a unified manner via a generic 1-D transport equation for horizontally averaged property. A water body interacts with underlying sediments. These sediments are represented by a set of vertical columns with heat, moisture and CH4 transport inside. The model is validated vs. a comprehensive observational data set gathered at Kuivajärvi Lake (southern Finland), demonstrating a fair agreement. The value of a key calibration constant, regulating the magnitude of methane production in sediments, corresponded well to that obtained from another two lakes. We demonstrated via surface seiche parameterization that the near-bottom turbulence induced by surface seiches is likely to significantly affect CH4 accumulation there. Furthermore, our results suggest that a gas transfer through thermocline under intense internal seiche motions is a bottleneck in quantifying greenhouse gas dynamics in dimictic lakes, which calls for further research.

The other mode of freshwater bodies impact on climate, is that through emissions of carbon dioxide (CO 2 ) and methane (CH 4 ) to the atmosphere (Tranvik et al., 2009). For instance, according to recent estimates (Bastviken et al., 2011), global methane flux from lakes offsets 25% of the esti-1 mated land carbon sink, implying that lakes are an important component of global carbon cycle and 20 climate system.
Concomitantly with growing awareness of lakes significance for current and future climate change, only few attempts have been made to develop lake models, incorporating thermodynamics, turbulence and biogeochemistry in order to simulate methane and carbon dioxide in natural water bodies (Stepanenko et al., 2011;Kessler et al., 2012;Tan et al., 2015). The ultimate goal of these devel-25 opments is to study the response of lakes and their greenhouse gas emissions under future climate change scenarios (Tan and Zhuang, 2015b) and through implementation of biogeochemical lake models into the Earth system models.
However, a number of problems arise concerning CH 4 and CO 2 modelling in lakes. First, variety of biogeochemical processes involved in production and transformations of CH 4 and CO 2 30 are not well understood to an extent where rigorous mathematical description could be developed.
For instance, methane production dependence on environmental factors has been tested in a bulk of studies (Borrel et al., 2011), however, to the best of our knowledge, only temperature dependence is quantified with high statistical confidence, e.g. (Yvon-Durocher et al., 2014). Moreover, even a widely-accepted statement that methane is produced exclusively in anaerobic environment 35 faces contradiction with some observational results (Damm et al., 2010), suggesting that there are CH 4 production mechanisms, not comprehended so far even at qualitative level. Second, lakes vary very much in climate, geological and biogeochemical environment resulting in enormous variability in greenhouse gas status (Juutinen et al., 2009). This situation is complicated by high vertical and sometimes horizontal variability of gas concentrations in a given lake (Schilder et al., 2013;Blees 40 et al., 2015). Third, when considering gas dynamics in lakes new physical processes become crucial, such as diffusion through the water surface (Donelan et al., 2002), vertical diffusion in metalimnion and hypolimnion, bubble interactions with sediments skeleton (Scandella et al., 2011), and others.
Many of these have not been addressed enough so far in both theoretical and experimental studies.
The obstacles described above hinder development of mathematical model on fundamental process-45 understanding basis. Therefore, any lake greenhouse gas model would contain a number of empirical constants to be calibrated on an extensive dataset (Tan et al., 2015), what is a usual practice in e.g. wetland methane models (Walter et al., 1996;Walter and Heimann, 2000;Wania et al., 2009;Melton et al., 2013).
This work aims at developing the lake model based on rigorous mathematical development feasi-50 ble in framework of 1D approach, applied for thermodynamic, hydrodynamic and biogeochemical prognostic variables in unified manner. We avoid using procedures for formal optimization (calibration) of the model parameters, rather focusing on qualitative behaviour of the model and its sensitivity to selected uncertain processes and constants. Vertical turbulent flux of dissolved gases through hypolimnion and metalimnion are of special concern in this work. The lake model, developed here 55 is based on LAKE model, that has been continuously advanced during last decade in Moscow State 2 Geosci. Model Dev. Discuss., doi:10.5194/gmd-2015-261, 2016 Manuscript under review for journal Geosci. Model Dev. lake. Then, Section 4 presents comparison of model results to observed data in a reference model run. In Section 5, we analyze results of reference experiment as well as of sensitivity experiments, elucidating significance of certain physical processes. Conclusions are summarized in Section 6.
2 The model overview LAKE model is a one-dimensional model solving horizontally averaged equations for heat, gases 70 and momentum transport for an enclosed water body. For taking into account heat and gases exchange with sloping bottom, the scheme for water temperature and gases concentrations is coupled to sediments columns originating at the bottom at different depths (see Section 2.5). Below we provide basics of 1D approach used and a general description for main groups of processes represented in the model. 75

The generic 1D equation and vertical coordinate
We commence the description of LAKE model with derivation of a generic 1D lake modelling framework, implemented in current version of the model in respect to all prognostic variables. We confine ourselves to a concise summary of that derivation, while the interested reader will find a rigorous mathematical development in Appendix A. 80 We start with the generic Reynolds-averaged advection-diffusion equation for the quantity f , that might be one of horizontal velocity components, temperature, turbulent kinetic energy (TKE), TKE dissipation or gas concentration (hereafter using summation over repeated indices): assuming mass conservation equation for incompressible fluid: where u i is the velocity component along x i Cartesian axis (x 3 = z being an axis pointing along gravity and originating at a lake surface, x 1 = x, x 2 = y -horizontal coordinates, u 1 = u, u 2 = 3 Geosci. Model Dev. Discuss., doi:10.5194/gmd-2015-261, 2016 Manuscript under review for journal Geosci. Model Dev. with A(z) denoting the area of horizontal cross-section of a lake at depth z. After applying this operator to (1) and making use of appropriate simplifications (Appendix A) we get: where we have decomposed the total vertical flux F z = F 3 into turbulent flux, F tz , and a nonturbulent flux, F nz , F z = F tz +F nz ; u h = (u 1 , u 2 ), n being an outer normal vector to the boundary Γ A(z) of cross-section A(z), k f -turbulent diffusion/dissipation for variable f , and a subscript b indicating a variable's value at the sloping bottom. The vertical fluxes of quantity f at lake's margins, 100 F tz,b , F nz,b , we will thereafter call marginal fluxes for brevity (marginal heat flux, marginal gas flux, marginal friction, etc. For the horizontally mean turbulent flux we applied a first-order closure, The non-turbulent fluxes enter equations for temperature (shortwave radiation flux) and for gases' concentrations (bubble flux).
when applying the lake model for the lake under study, we will use Wedderburn number time series to check the validity of dropping out the "vertical circulation term" 1 .
Equation (4) is a generalization of equations that include lake shape effects encountered in many 1D models designed for lakes (Stefan and Fang, 1994;Goudsmit, 2002;Jönhk et al., 2008;Tan and Zhuang, 2015a) as well as for reservoirs (Zinoviev, 2014). In all 1D lake models we are aware of, 130 term IV does not include shortwave radiation flux in temperature equation and misses bubble flux of gases in equations for dissolved CH 4 and CO 2 .
The form of equation (4) written using geometric vertical coordinate z is not convenient for the case of significant rate of water level change. In order to tackle this case, a normalized vertical coordinate, ξ = z/h(t), where h is the maximal depth of a lake, has been introduced into equations 135 of the model. Furthermore, the movement of z-axis origin when the surface level changes, strictly speaking, results in an additional term to a generic equation (4). Above leads to the following final form of (4): with B s signifying precipitation minus evaporation, i.e. the rate of z-axis origin motion, positive 140 upwards. Although, it is the form (5) that is implemented in the LAKE model, it differs from (4) by metric terms only, so that for the sake of simplicity in subsequent flow we will refer to equation (4).
Moreover, in this work we will keep lake depth h constant, realistic for the lake under study.
1 Other possible mechanisms for basin-scale circulations include density currents along sloping bottom (Chubarenko, 2010;Kirillin et al., 2015) during transitional seasons and ice period.

Lake thermodynamics
The water temperature in the model is driven by equation (4) with substitution f → T , where 145 c = c w ρ w0 , c w -water specific heat, ρ w0 -reference water density, R f = 0 (no heat sources in the water besides radiation heating), F nz (z) = F nz,b (z) = S rad -shortwave radiation flux, positive downwards. The latter equality means that we assumed shortwave radiation flux to be horizontally homogeneous at all depths. This is commonly used approximation as getting data of the spatial distribution of turbidity in a lake requires special measurements. Heat conductance is a sum of molecular 150 and turbulent coefficients, Shortwave radiation flux, S, is treated as consisting of near infrared fraction and the rest energy (mostly visible radiation). Near infrared part is consumed completely at the surface, whereas the visible fraction is partially reflected according to water albedo, and its remainder is attenuated with 155 depth according to widely-used Beer-Lambert law with extinction coefficient specific for a lake under study (see Section 3.2).
To solve (4) for temperature one needs to specify top and bottom boundary conditions as well as a method for calculation of marginal heat flux, F tz,b (z), at each depth z. The top boundary condition is a well-established heat balance equation, involving net radiation and a scheme for turbulent heat 160 fluxes in a surface atmospheric layer based on Monin-Obukhov similarity theory (Paulson, 1970;Businger et al., 1971;Beljaars and Holtslag, 1991). The way of coupling the water column to bottom sediments through lower boundary condition and marginal heat flux is less straightforward. When the heat transfer in bottom sediments is solved by diffusion-lake equation, there are two options for imposing boundary conditions at the "water-sediments" interface:

165
-Continuity of both heat flux and temperature at the interface; -Continuity of heat flux across the interface and a method for heat flux calculation, relating it to in-water temperature gradient, e.g. through logarithmic profile formulae.
We found that the first option provides reasonable results for temperature and especially for gas concentrations (see below in the paper), whereas the second one needs calibration of parameters 170 entering the flux-gradient relationship in the bottom boundary layer. The marginal heat flux is calculated using the same temperature and flux continuity condition, that is facilitated by the solution of vertical heat transfer in sediments below sloping bottom (see details in Section 2.5).

Lake hydrodynamics
Applying the form (4) to horizontal momentum equations is straightforward with F nz = 0, c = 1 and R f representing Coriolis force and horizontal pressure gradient. The Coriolis force has to be included in the momentum equations for lakes with horizontal size exceeding internal Rossby deformation radius (Patterson et al., 1984), that we will check below when validating the model for 180 the lake under study.
The term F tz,b (z)A −1 dA/dz has the sense of marginal friction for the case of momentum equations. This term can be parameterized as quadratic in velocity with tunable proportionality coefficient (Jöhnk, 2001). Instead, we apply logarithmic layer friction with effective bottom roughness length, z 0b,ef f . The characteristic "effective" in respect to z 0b,ef f accounts for the fact that while calculat-185 ing bottom friction we use horizontally averaged velocity components u, v instead of the velocity components' values in the logarithmic layer adjacent to bottom. As there are no theoretical hints how z 0b,ef f relates to the "true" bottom roughness, z 0b , it may be used as tunable parameter. However, our modelling results show, that choosing z 0b,ef f of the order of z 0b expected at the bottom eventually provides reasonable results in terms of vertical mixing of water properties.

190
The more interesting story comes with parameterization of horizontal pressure gradient. We represent it at any depth z as (h s -lake surface deviation from horizontal), implying that we have used hydrostatic equation with constant density, ρ w0 . This is a barotropic approximation since we neglected buoyancy in the hydro-195 static equilibrium 2 . It is the simplest way to account for horizontal pressure differences still being capable to induce significant mixing below the thermocline (see below in the paper).
To estimate terms (6) in 1D model we modify the scheme proposed originally by U.Svensson (Svensson, 1978;Goudsmit, 2002). Fig. 1 where L x , L y are the sizes of horizontal water body intersection in x and y directions, respectively, 205 subscript "0" denotes values at the lake surface. For simplicity, in the model we approximate the lake's horizontal cross-section, A(z) as an ellipse, so that L x stands for twice of major semiaxis and L y -for twice of minor semiaxis, or vice versa. Equations (9)-(10) express the change of surface level of four lake's sections (h s,x1 being the mean of h s over the "left" section of a lake, x < x c , h s,x2 -the same for the "right" lake section x > x c , h s,y1 -for y < y c , h s,y2 -for y > y c , (x c , y c ) 210 standing for the lake center) due to volume discharge through two vertical planes, x = x c and y = y c , ( Fig. 1), neglecting inflows and outflows. The multiplier π 2 /4 in (7)-(8) arises instead of a "natural" choice of 2 in order the solution of the model equations for specific case of rectangular channel be matching the period of the 1-st surface seiche mode, i.e. the Merian formula (Merian, 1828) (see Appendix C for mathematical development). According to that, the output of LAKE model for the 215 case of 1D flow developing along a non-rotating channel after initial disturbance of lake surface demonstrates oscillations with a period very close to that predicted by Merian formula (not shown here). In the following, we will refer to 7-10 as either "surface seiche parameterization"or "dynamic pressure gradient parameterization".
Boundary conditions for momentum equations are momentum flux from the atmosphere, calcu-220 lated according to air surface layer bulk formulae (Paulson, 1970;Businger et al., 1971;Beljaars and Holtslag, 1991), and friction at the deepest part of bottom following quadratic dependence on velocity with Chézy. Momentum flux accelerating currents is parameterized as a fraction of total momentum flux from the atmosphere (Stepanenko et al., 2014), because in conditions of limited fetch (small lakes) a part of total momentum flux is consumed by wave development. Partitioning mo-225 mentum flux between waves and in-water currents significantly reduces shear-driven vertical mixing during summertime.

Turbulence closure
The turbulence closure is a k − model with Canuto stability functions (Canuto et al., 2001). Nonturbulent flux, F nz in (4) is put to zero, because this model does not include any fluxes of k and 230 besides advection and turbulent transport. We neglect also advection of TKE and dissipation rate  (Burchard and Petersen, 1999). Sinks and sources of TKE and dissipation rate, i.e. buoyancy term and shear production, hidden in R f of (4), are approximated using only vertical derivatives of horizontally averaged temperature, salinity 240 and velocity components. For constants of k − model used in this study, see Appendix B.
In our study the turbulence closure briefly described above will be referred to as "standard k − model". Pertinent to objectives of the study, we will also use extensions of standard k − model to account for specific mixing mechanisms in the thermocline, namely, gravity waves (Mellor, 1989) and internal seiches (Goudsmit, 2002).

Heat and moisture processes in sediments
Snow and ice modulae are not used in this study. Processes in sediments are treated inside a set of 3D figures, that all have the same vertical dimension, h sed , and the horizontal intersections of which are confined by sequential isobaths (Fig. 2). In each such a column all properties of sediments are assumed to be horizontally homogeneous, so that only the vertical transport of heat and other 250 quantities applies. Each column of sediments exchanges heat and methane with the horizontal water layer bounded from below and above by respective isobath levels according to continuity of flux and a quantity considered (temperature, methane concentration, see Section 2.2).
The heat processes in the model include vertical transport and phase transition between water and ice. The vertical transport in sediments is described according to (Côté and Konrad, 2005).

255
Liquid water is transported via gravity and capillary-sorption forces (Stepanenko and Lykossov, 2005). The latter are represented by diffusion-like term. Bottom boundary condition for temperature is geothermal heat flux, usually set to zero. For moisture, saturation of sediments is used for the top boundary and zero flux is applied at the bottom.
2.6 Biogeochemistry and transport of CH 4 , CO 2 and O 2

Methane in sediments
Elaborate description of CH 4 model in sediments can be found in (Stepanenko et al., 2011), whereas here we provide a general overview and latter amendments to the model presented therein. This 265 model is applied in every column of sediments under a lake.
In each column of sediments (Section 2.5) methane transport is considered to be vertical only. The governing equation for bulk methane concentration, C CH4 , reads: where k CH4,s designates molecular diffusivity of methane, P CH4,s -production rate, E CH4,s -ebul-270 lition rate, O CH4,s -aerobic oxidation rate (anaerobic oxidation is omitted), and z s denotes a vertical coordinate originating at the column top. Vegetation uptake of methane by roots and aerenchyma transport are neglected in this study. Methane production rate is confined to the upper part of a sediments column and controlled by temperature by exponential dependence:

275
with P 0 being a calibrated constant, reflecting the amount and quality of organic material in sediments in respect to methane production, α new = 3 m −1 -a constant controlling the decrease of CH 4 production with depth, H is a step (Heaviside) function, q 10 = 2.3 (Liikanen et al., 2002) -temperature dependency constant, T mp -melting point temperature, α O2,inhib -a constant describing the rate of inhibition of methane production with rise of bulk oxygen concentration in sediments, C O2,s .

280
The latter constant is set as α O2,inhib = 316.8 m 3 /mol to ensure 100 times inhibition of methane production at oxygen content of 10 ppm, implying almost complete suppression of methanogenic Archaea activity under this concentration (Borrel et al., 2011). Parameterization (12) traces back to (Walter et al., 1996), the last multiplier added in this study. Deep methane production from old organics (Stepanenko et al., 2011) at the bottom of talik is not included as the lake under study is not 285 a thermokarst one.
Ebullition rate, E CH4,s becomes non-zero when bulk methane concentration exceeds a critical value, defined by hydrostatic load of water column and sediments layer above at a given depth, z s , as well as by nitrogen concentration at the sediments top (Stepanenko et al., 2011;Walter et al., 1996). Retention of bubbles in a sediments skeleton (Scandella et al., 2011) is neglected, so that 290 methane ebullition flux at the sediments top of k-th column, F B,1,k , is calculated as 300 where V max,s = 1.11 * 10 −5 mol/(m 3 * s), K CH4,s = 9.5 * 10 −3 mol/m 3 , K O2,s = 2.1 * 10 −2 mol/m 3 are methane oxidation potential and two half-saturation constants, respectively (Lidstrom and Somers, 1984).
In order the above scheme for methane oxidation and methane production inhibition to be realistic, the top numerical layer in sediments is set to be of thickness typical for oxygenated layer in lake's 305 sediments, 1 cm (Huttunen et al., 2006).

Methane in water
Methane concentration in water evolves according to equation of the form (4) with the term I representing the input of CH 4 by inlets and its outflow by outlets (not taken into account in this study). Diffusion coefficient, k CH4,w , is set equal to heat conductivity (turbulent Lewis number

Oxygen and carbon dioxide in water
Oxygen concentration is simulated by equation (4) with term I neglected, assuming turbulent Lewis number to be 1, and marginal diffusive flux treated as sedimentary oxygen demand (SOD). Other sinks of oxygen are biochemical oxygen demand (BOD, excluding respiration), respiration and 320 methane oxidation. Methane oxidation bacteria consume oxygen according to widely-accepted stoichiometric relation providing the rates of oxygen consumption and CO 2 production given the rate of methane loss (Section 2.6.2). The only process producing oxygen in a water column is photosynthesis. For bio-325 chemical oxygen demand, respiration and photosynthesis we use parameterizations from (Stefan and Fang, 1994). These parameterizations assume the rates of biogeochemical processes to depend exponentially on temperature and be proportional to chlorophyll-a concentration. Photosynthesis is additionally limited by photosynthetic active radiation. In our simulations, we kept the original empirical constant values from (Stefan and Fang, 1994). For more details an interested reader may refer 330 to the original paper.
As to sedimentary oxygen demand, we adopted the formulation from (Walker and Snodgrass, 1986), as it involves explicitly the near-bottom oxygen concentration (via diffusive term), in contrast to that from (Stefan and Fang, 1994), where SOD continues to be non-zero even when oxygen content in water is nil.

335
Carbon dioxide in water is calculated by the same type of prognostic equation as for other gases.
The only sink of CO 2 in the water column is photosynthesis, whereas its production in the model is provided by SOD, BOD, respiration and methane oxidation. As the rates of these processes in terms of O 2 and CH 4 are quantified above, the respective income or loss of CO 2 is immediately provided by (15) and the following stoichiometric equality: 2.6.4 Diffusive gas flux at the water-air interface The top boundary condition (at the lake-atmosphere interface) for concentration of any dissolved gas has the form: where C w is C CH4,w , C O2,w or C CO2,w , k s -dissolved gas diffusion coefficient and F Cw is the diffusive flux of a gas into the atmosphere, positive upwards. This flux is calculated according to the widely used parameterization: with C ae being the concentration of the gas in water equilibrated with the atmospheric concentration following Henry law and k ge , m/s, denoting the gas exchange coefficient, the so-called "piston velocity". The latter is written as: with the Scmidt number Sc(T ) having individual values for different gases and being temperature-355 dependent. The k 600 coefficient has been a subject of numerous studies, and a number concepts have been proposed to quantify it (Donelan et al., 2002). We adopt two widespread options for k 600 : (i) empirical dependence on wind speed and (ii) surface renewal model.
The dependency on wind velocity takes the form (Cole and Caraco, 1998): Here, u a,10 stands for the wind speed vector at 10 m above the water surface, C k600,1 = 5.75 * 10 −6 m/s, C k600,2 = 5.97 * 10 −6 (m/s) 1−nk 600 are empirical constants. The simple empirical equation (21) "integrates" the effects of wind speed on a number of processes such as turbulence in adjacent layers of water and air, wave development and breaking, cool skin dynamics, and therefore is likely to be not enough sophisticated to express adequately a wide variety of conditions met on 365 real lakes. Therefore, we also included surface renewal model (MacIntyre et al., 2010;Heiskanen et al., 2014), that in terms of k 600 states: where ν m designates molecular viscosity of water, and C 1,SR = 0.5 is an empirical parameter. As TKE dissipation rate is available directly from k − closure, we do not use any special parameteri-370 zation for | z=0 as proposed in other works (e.g., (MacIntyre et al., 2010)).

Single bubble model
The bubble model used in LAKE closely follows that described in (McGinnis et al., 2006). Consider the evolution of a bubble rising from the lake bottom, and consisting of a mixture of gases.

375
The quantity of each i-th gas in the bubble, M i , mol, changes due to its dissolution into the water according to equation: where r b , m, is the bubble radius, H i -the Henry "constant" dependent on temperature T , K, P i , Pa, the partial pressure of i-th gas, C i , mol/m 3 , is the concentration of a gas in water, K i is the at the bottom and pointing opposite to gravity, n g is the number of gases in a mixture.
Five gases are considered in a bubble: methane, carbon dioxide, oxygen, nitrogen and argon.
Water vapour constitutes minor contribution to a bubble pressure, and therefore neglected. Indeed, the saturated vapour pressure at 20 • C is 23.4 hPA that is ≈ 2% of atmospheric pressure. This is 385 the upper estimate for water vapour pressure contribution in bubbles, as the pressure increases with depth, and saturation vapour pressure -decreases, due to water temperature decrease.
The temperature in the bubble is assumed to be equal to that of environmental limnetic water at the depth of current bubble location, Z. It means that the heat exchange between the rising bubble and water is expected to be intensive enough to dominate over the adiabatic cooling of the bubble. In 390 practical terms, this frees us from solving additional equation for bubble temperature. The temperature dependency of Henry constants for flat solution surface is taken from (Sander, 1999). The effect of gas-water interface curvature on equilibrium gas pressure is omitted in this model because when using Thomson (Kelvin) formula it turns out to be negligible for typical bubble radii in oceans and lakes (≥ 1 mm). Exchange coefficient, K i , is dependent on molecular diffusivity in water, bubble 395 radius and its velocity according to empirical formulae from (Zheng and Yapa, 2002). The bubble velocity is determined assuming equilibrium between buoyancy force and environment resistance given by quadratic law for small radii (r b < 1.3 mm) and taking into account bubble surface oscillations for larger sizes (Jamialahmadi et al., 1994).
For each component of gas mixture we apply an ideal gas law because under the typical pressures 400 at moderate water depths (at least dozens of meters) Van der Waals forces are small: where R is the universal gas constant. The surface tension pressure is small for the bubbles with radii typical in lacustrine environment, and is neglected in (24). Then, when equating the gas mixture pressure ng i=1 P i to hydrostatic pressure at a given depth, p a + ρ w0 g(h bot − Z) (p a the atmospheric 405 pressure, h bot is a lake depth in a point, where the bubble is released) and using (24) one yields: where M 0 -the total gas quantity in the bubble (mols). According to (26), the bubble initialization is provided by initial bubble radius, r b0 , and molar fractions of mixture components α i . In this study, we chose r b0 = 2 * 10 −3 m and the initial bubble gas composition to be 100% of CH 4 .
The bubble model described above is numerically solved by Euler explicit scheme.

Bubble flux of gases
In equation (4) applied for methane, oxygen and carbon dioxide the non-turbulent flux (term III) consists of bubble flux only. Bubble flux also contributes to term IV therein. This section explains how these terms are evaluated using single bubble model, described above (Section 2.7.1).
We consider an idealized situation when all bubbles rising from all columns of sediments have 420 the same initial radius r b0 at the bottom and identical gas composition. Given that in reality there is always a distribution of bubbles in their size, parameter r b0 may be treated as an average (in appropriate sense) radius over this distribution. For bubbles rising from a given sediments column, equations (23)-(25) imply that their radius and composition will be the same at any level over this column.

425
Now, at any depth z we can construct a horizontal average of vertical bubble flux of i-th gas, where index k is the index of a sediments column, n s being a total number of columns, and A k (z) is an area of projection onto A(z) of the part of the top facet of k-th column residing below depth z (e.g., for columns with tops above z, A k (z) = 0; for columns of sediments with top facets completely below z, A k (z) = A s,k , A s,k standing for the area 430 of top facet of k-th column). When the mean flux is calculated, it may be used in the term III of (4) Here F B,i is defined as positive upwards leading to "+" sign.
To get the averaged flux F B,i as described above, individual bubble fluxes F B,i,k are calculated from each sediments column as Here, we introduced bubble number density n b,k , m −3 , and k is a sediments column index, as before. All bubbles that are released from a given sediments column's surface completely dissolve simultaneously at some depth or evade to the atmosphere. Furthermore, it is known that bubbles with diameter ≈ 1 cm are unstable and split up (Yamamoto et al., 2009;McGinnis et al., 2006). Hence, 440 in the model it is assumed that a bubble with r b ≥ 0.5 cm splits into two. In the depth interval between two subsequent bubble collapses the bubble flux (that is the number of bubbles crossing the horizontal surface of 1 m 2 per second) is constant, and at the depth of division it doubles. Taking this into account, one may rewrite (28) as follows

Numerical aspects
The principal requirements for the numerical scheme of the diffusion-type model with nonlinear sources described above are integral conservation of prognostic variables and stability.

455
Integral conservation is achieved by employing second-order centered differences in space for all equations in water and sediments. The coupling of sediments columns to water body is also implemented ensuring continuity of heat and methane flux across the sediments-water interface.
Equations of k − closure are discretized in a way where TKE input by shear production and buoyancy in TKE equation equals to dissipation and potential energy source/sink in momentum and 460 temperature/salinity equations, respectively (Burchard, 2002) (salinity is set to zero in this study).
Time-marching scheme is a Crank-Nicolson one (Crank and Nicolson, 1996) that allows for increased time steps, ∆t ≈ 10 min for vertical grid spacing of ≈ 1 m in water, if not using surface seiche parameterization. The time step is limited due to high non-linearity of k − closure. However, the strongest constraint for time step arises when horizontal pressure gradients are calculated 465 via mass conservation (7) The parameters of baseline experiment are given in Table 1. Maximal lake depth was set to 12.5 m to ease comparison with measurements, as this is the local depth below observational mast. There is no information on the lake sediments characteristics of Kuivajärvi Lake, however, silt loam should be close in soil particle size to typical lake sediments. Sediments depth (10 m) is chosen to be of enough extent for temperature fluctuations not reaching its lower boundary. To get A(z), we linearly 510 interpolated the morphometric data given in Table 2.
Boundary conditions were set as follows. At the sediments bottom, zero heat and moisture fluxes were imposed. At the water surface, heat balance equation is applied, where downward radiation fluxes were measured at the mast, surface longwave radiation calculated via Stefan-Boltzmann law, and sensible and latent heat fluxes -using Monin-Obukhov similarity functions (Table 1). In total, 515 seven meteorological variables were supplied to the model, at 30 min interval, all measured at the lake: wind speed and direction, temperature, humidity, longwave and shortwave radiation, atmospheric pressure. For analysis of these time series, we refer to .
Initial conditions for the model are the profiles of all prognostic variables at initial instant. Water temperature, oxygen, carbon dioxide and methane vertical distributions were specified from mea-520 surements at 00:00 5 May 2013. Salinity was set to zero, two horizontal velocity components were initialized with small values. In sediments columns, temperature was set to 4 • C, and water content -to slightly undersaturated values.
Only two constants in the model were calibrated. The first one is P 0 in equation (12), controlling the magnitude of methane production in sediments, representing quantity and "quality" of organics 525 in sediments as a substrate for methanogenic activity. However, we found that this constant is not enough to regulate methane concentration in the lake mixed layer (see the rationale in Section 5.3). A half-saturation constant in CH 4 oxidation reaction rate, K CH4,w , was found to be crucial parameter in this respect, effectively changing mean levels of mixed-layer methane concentration.
The sensitivity experiments were set with the same configuration, as the baseline experiment, with  Initial bubble radius, rb0 2 * 10 −3 m

Physical parameterizations
Surface flux scheme (Paulson, 1970;Businger et al., 1971;Beljaars and Holtslag, 1991) Equation of state (Hostetler and Bartlein, 1990) Turbulence closure standard k − with Canuto stability functions (Canuto et al., 2001) minimal diffusivity in the thermocline increased (MD) In the following sections we will describe and discuss main results of the baseline experiment and sensitivity experiments in terms of physical and biogeochemical variables.

Temperature and turbulent quantities
In this section we will consider temperature stratification and turbulent structure of the lake vertical column, that are prerequisites for correct simulation of biogeochemical processes. The surface momentum and energy fluxes will not be covered as they were discussed for this lake involving LAKE 545 model results on these variables in .
Evolution of temperature distribution in the lake is presented at Fig.5a (model) and  The model satisfactorily reproduces the observed seasonal temperature pattern in Kuivajärvi Lake.
The root mean square error (RMSE) for surface temperature is 1.54 • C and the difference of means is 0.61 • C (the average of modeled surface temperature slightly exceeds that of observed). However, a closer look into Fig. 5b reveals high-frequency fluctuations of observed temperature in the depth range of the thermocline which are not reproduced by the model (Fig. 5a). These oscillations are of

570
On the other side, Mellor gravity waves parameterization (Mellor, 1989) (GV+ experiment) adds considerable TKE to the profile of baseline experiment, especially in metalimnion and hypolimnion in conjunction with surface seiches taken into account.

Oxygen
Hereafter, for gases dissolved in water, we use concentrations per unit volume of water.

575
Dissolved oxygen evolution in Kuivajärvi Lake is presented at Fig.7a (model, reference experiment) and Fig.7b (observations). Here and in the subsequent plots for CO 2 and CH 4 , we will confine our analysis to June-October period, as during May the modeled gas concentration undergo adjust- ments column locating approximately below the EC footprint that should be used (for EC footprint at Kuivajärvi Lake, see ). In our case, it is the deepest sediments column, where the time average CH 4 ebullition flux reaching the surface constituted 0.006 µmol/(m 2 s) (8 mg/(m 2 day)). Thus, the mean total methane flux in the model exceeded the observed one by an order of magnitude, still remaining low compared to that at many other lakes (Juutinen et al., 2009). lake's layered structure, is assessed via Wedderburn number (W ) (Shintani et al., 2010) or Lake number (Imberger and Patterson, 1989), while the significance of Coriolis force can be quantified by comparing Rossby deformation radius (L R ) to a lake size (Patterson et al., 1984). The two parameters, W and L R , are plotted in Figs. 10 and 11, respectively. We see that during June-August W generally fluctuates around 50 implying that thermocline vertical displacement by wind forcing 635 is about ∼ 100 times less than the mixed-layer depth. However, in May, end of October and on short periods of several days during summer W approaches unity making possible upwelling at lake's margins, similar to what was reported for 2011 in (Heiskanen et al., 2014). At least two of these episodes, namely, these in mid-June and end of July, are concomitant to mixed-layer cooling (Fig.   5a), weakening the lake stratification.  (Glazunov and Stepanenko, 2015) or laboratory experiments (Markfort et al., 2013)). However, it turns out in practice of 1D lake model applications to such lakes, that not only Monin-Obukhov is the most physically-based 655 option in these models for obtaining surface heat fluxes so far, but it is that still delivering acceptable accuracy for calculated surface temperature at seasonal timescales (Stepanenko et al., 2014;Heiskanen et al., 2015). As to momentum flux, it has been shown to be a crucial parameter regulating the rate of mixed-layer deepening during summer (see e.g. (Stepanenko et al., 2014)). This resulted in a widespread modelling practice where the drag coefficient, defining momentum flux at a given wind-660 speed, has become a tunable parameter in k − -based lake models. In our simulations, we do not calibrate surface drag coefficient, but include a simple parameterization of momentum flux partitioning between waves and currents (Stepanenko et al., 2014)  and ω ≈ 4.5 * 10 −4 s −1 (T seiche ≈ 3.9 hour). The harmonic of T seiche ≈ 20.5 hour contains much more energy than that of T seiche ≈ 3.9 hour. In order to interpret these spectra we use the method for seiche period calculation proposed by (Münnich et al., 1992).
Starting from two-dimensional linearized incompressible Boussinesq equations and seeking the 680 solution for vertical velocity, w, in a wave-like form: with rigid lid condition w| z=0,h = 0 leads to an ordinary differential equation for the amplitude, W :

685
which is a Sturm-Liouville problem for frequencies, ω and corresponding eigenfunctions, W . We solved it by shooting method with squared Brunt-Väisälä frequency N 2 taken from the mean temperature profile measured in July and h = 12.5 m (a depth of lake in the point of measurements).
Considering 1-st horizontal mode, k = π/L x0 , we've got T seiche,1 = 6.5 hour with W having a form of the 1-st vertical mode, usually denoted as V1H1 (one maximum of W between z = 0 and 690 z = h) and T seiche,2 = 21.2 hour for the 2-d vertical mode (one maximum of W and one minimum, V2H1). These frequencies correspond to those of maxima at the temperature spectrum (Fig.12).
The discrepancy between measured and calculated frequencies, that especially noticeable for V1H1 mode (3.9 hour vs. 6.5 hour, respectively), is expectable since the linear analysis described above neglects morphometry of the lake's bed (Fricker and Nepf, 2000), effects of Coriolis force and the 695 complex temporal behaviour of the actual wind forcing.
The prominence of V2H1 mode in the temperature spectrum is what have been found for an Alpine lake by M.Münnich as well (Münnich et al., 1992). A plausible explanation for that is the resonance between V2H1 seiche and the wind speed, both having close to diurnal periodicity (Mortimer, 1953).
Thus, the main conclusion of this section is a presence of significant internal seiches in Kuivajärvi 700 Lake that may be responsible for additional mixing in the thermocline either in the interior of the lake or at its margins. This will be discussed in the following section (Section 5.2).

Turbulent quantities
In this section we will focus on turbulence characteristics in the thermocline and hypolimnion as they are factors for vertical transport of gases originating at a lake bottom. Moreover, the presence 705 of seiches in the lake suggests additional mixing mechanisms to exist in the thermocline, such as production of TKE by near-bottom shear (Goudsmit, 2002) and breaking of internal waves at the sloping bed (MacIntyre et al., 2009;Boegman et al., 2005).

TKE production terms
The vertical distribution of TKE shown at Fig.6 is formed as a result of approximate balance between 710 terms in right-hand side of TKE equation (B1). Mean vertical distribution of TKE production by shear, S, by buoyancy, B, and by seiches, S seiche (only when Goudsmit parameterization is used, "IS+" experiment) in July is shown in Fig.13.
First, we see that mean buoyancy production is positive in the top half of mixed layer (∼ 10 −9 ÷ 10 −7 m 2 s −3 ), indicating that nocturnal buoyancy production of TKE in this region overrides the 715 daytime sink. It is several times (up to an order of magnitude) less than the shear production, however, exceeds 3-5 orders of magnitude generation of TKE by seiches. Different experiments show almost identical profiles of B. It is because both Goudsmit and Mellor parameterizations include dependence on N providing zero contribution to TKE and other turbulent quantities at N = 0, and N ≈ 0 in the mixed layer. Below, buoyancy production becomes negative due to stable stratification.

720
Vertical shear production is the largest contributor to TKE throughout a lake profile excepting thermocline, at ∼ 7 m depth, where it attains its minimum and becomes less than S seiche . This minimum corresponds to TKE minimum (Fig.6) and a minimum of eddy viscosity, ν, approaching minimal value, ν min . As in the model we do not use any "background diffusivity/viscosity/conductivity", the minimum value of ν and ν T is set to a very small number, ν min = ν T,min = 10 −8 m 2 /s (cf. molec-725 ular viscosity at 10 • C, ν m = 1.307 * 10 −6 m 2 /s and heat diffusivity, ν T,m = 1.41 * 10 −7 m 2 /s).
Hence, S = ν[(∂u/∂z) 2 + (∂v/∂z) 2 ] reaches negligible values, as ν = ν min . Below thermocline, there is drastic difference in S between experiments where dynamic barotropic pressure gradient was taken into account (baseline experiment), and those without surface seiches -labelled by "SS-" on Fig.13. The reason is that due to water surface inclination, currents are generated in hypolimnion, shear production takes place in the experiment with both Mellor parameterization and dynamic pressure gradient included ("GV+"). The value of S is especially increased in the thermocline, because N 2 reaches maximum there, and it contributes to corresponding additional shear proportionally.
TKE also achieves maximal values for this experiment at all depths (Fig.6). Still, heat conductance 735 and diffusivity in the metalimnion are close to molecular values even in this case.
Additional shear production due to seiches attains maximum in the thermocline with minima in the epilimnion and hypolimnion. This is again due to proportionality S seiche ∝ N 2 . The contribution of S seiche to TKE production remains minor compared to shear everywhere, excepting a small region in the thermocline, where TKE generation by vertical shear plunges to minimum, as discussed above.

740
The strong effect of surface seiches on under-thermocline turbulence obtained in our study is yet to be verified with more complicated models (i.e. 3D Reynolds-Averaged Navier-Stokes or Large Eddy Simulation) and extensive turbulence measurements for Kuivajärvi Lake. Indeed, surface seiches are barotropic motions not taking into account density stratification in a lake. As a consequence, their period of ∼ 1 min for Kuivajärvi Lake is orders of magnitude less than that of V1H1 mode (6.5 h) 745 and higher modes obtained from eigenvalue problem for continuous stratification (Section 5.1). Taking into account internal seiches in the model would change drastically frequencies of near-bottom current oscillations compared to surface seiches and thereby the hypolimnetic shear production of TKE. However, so far, to the best of our knowledge, internal seiche parameterization producing extra mixing in the hypolimnion has not been developed, as the poineering attempt by (Goudsmit,750 2002) introduced S seiche ∝ N 2 , negligible in hypolimnion. Envisaging implementation of internal motions in the model for our future work we, however, note that introducing surface seiches allowed to generate TKE below thermocline qualitatively consistent with a bulk of observational data (Wüest et al., 2000;Wüest and Lorke, 2003), demonstrating that summer stratification in dimictic lakes is comprised of two turbulent layers disconnected by quasi-laminar thermocline.

770
On the other side, when k− model is supplemented by Goudsmit internal seiche parameterization (Goudsmit, 2002), i.e. when the shear production is modified as S * = S + S seiche , S = νM 2 , an expression for stationary Richardson number may be derived as well (see Appendix D): (u a , v a ) standing for wind vector in the surface layer, ν 0 -eddy viscosity at stationary turbulence 775 regime, C s -constant for a given lake including empirical parameters and lake morphometry characteristics. As there are no unique values of k, and ν 0 resulting from uniformity and stationarity conditions, we assume a small value ν 0 ≈ ν m leading to an upper estimate, Ri st = 0.30. Larger values of ν 0 , according to (36), would decrease Ri st .
Estimates provided above suggest that Ri st in k − model still remains under unity, when gravity 780 waves and internal seiche parameterizations are included. Thus, they cannot generate significant turbulence in the thermocline of Kuivajärvi Lake, where Ri >> 1. Indeed, in all experiments minimal eddy diffusivity in the thermocline was close to minimal possible one set in the code, 10 −8 m 2 /s, implying only molecular diffusion to perform vertical transport. Still, we envisage a possibility of mixing mechanisms rising total diffusivity above molecular levels in the metalimnion, given empir-785 ical evidences (e.g. (Saggio and Imberger, 2001)) and the fact that Mellor and Goudsmit parameterizations have not been tested thoroughly vs. extensive measurement data and/or LES and DNS simulation so far. Therefore, we conducted a sensitivity test on the influence of artificially increased diffusivity in the thermocline on gas concentrations ("MD" experiment, see next Section, 5.3).

Oxygen, methane and carbon dioxide 790
As we see at Figs 7a and 7b oxygen concentration is high in beginning of June not only in the mixed layer (8 − 9 mg/l), where it is produced by photosynthesis, but beneath the mixed layer as well (5 − 7 mg/l). This is due to maximal oxygen concentrations throughout a water column during the spring overturn in beginning of May. Afterwards, oxygen remains high in the mixed layer while it decreases to almost zero values in hypolimnion by August. Unforunately, so far, there is no observational data for Lake Kuivajärvi (e.g. turbulence measurements or any sediments data), facilitating to discern, whether it is enhanced turbulence below 820 thermocline and/or nearly homogeneous SOD distribution with depth that makes measured oxygen profiles to be much more even than these in the model.
Consistently, the same questions arise considering carbon dioxide distribution (Figs 8a and 8b).
Neglecting the spot of low CO 2 hypolimnetic amount in the measured pattern around mid-August, that might be due to measurement errors, we see larger uniformity in the measured vertical distribu-825 tion than in that calculated. Bottom concentration rises much faster in the model up to ≈ 16 mg/l by mid-August, whereas in observed field this level is attained by mid-September only. This fast bottom accumulation of calculated CO 2 corresponds to fast decrease of O 2 (Fig. 7a). This corroborates our suggestion above that either vertically even SOD or vertical mixing (or both) are misrepresented in the model under thermocline, as these processes affect CO 2 and O 2 in the way to homogenize 830 hypolimnetic profile.
We also note that abrupt increase of deep CO 2 concentration that took place according to measurements in September in the depth interval 8-12 m is absent in the model. We argue that this rise is unlikely to be caused by local aerobic decomposition of organic matter, as the oxygen is depleted 28 Geosci. Model Dev. Discuss., doi:10.5194/gmd-2015-261, 2016 Manuscript under review for journal Geosci. Model Dev.
near bottom by this time, < 1 mg/l (Fig. 7b), and this amount is far from enough to contribute 835 to CO 2 jump by 5 − 7 mg/l, given stoichiometric ratio or corresponding reactions O 2 : CO 2 ∼ 1.
Hence, we suppose this CO 2 is advected to the point of measurements from catchment. Moreover, this early autumnal sharp increase of carbon dioxide is likely to be a peculiarity of 2013; at least, in 2011 and 2012 rising CO 2 hypolimnetic concentration was much more smooth (Miettinen et al., 2015).

840
Kuivajärvi Lake is a significant source of CO 2 (Heiskanen et al., 2014;Miettinen et al., 2015;Mammarella et al., 2015), and significant underestimation of surface concentration by the model (0.39 mg/l vs. 2.80 mg/l measured) is a serious drawback of the model setup. As carbon dioxide in the mixed layer is affected by a large number of processes (BOD, SOD, respiration, photosynthesis, diffusion to the atmosphere), it seems for us difficult to disentangle this problem on a solid 845 physical/biogeochemical basis in this study, and it should be a part of separate research.
As stated above, only two constants were calibrated in the model, i.e. P 0 and K CH4,w , that are responsible for magnitude of methane production in sediments and methane oxidation in water, respectively. The value P 0 = 3 * 10 −8 mol/(m 3 * s) chosen occurred to be very close to the value obtained for thermokarst Shuchi Lake in North-Western Siberian (P 0 = 2.55 * 10 −8 mol/(m 3 * s), 850 see (Stepanenko et al., 2011)). We note that it is not straightforward to compare these values, because the model version used in Shuchi Lake study, lacked such important features as taking into account bottom morphometry, bubbles dissolution in water, all biogeochemical processes involving O 2 and CO 2 but methane oxidation. Nevertheless, the same order of magnitude of P 0 for two lakes of different genetic types with ecosystems functioning under drastically different climate conditions, 855 argues for robustness of our model formulation. A half-saturation constant for methane, K CH4,w = 3.75 * 10 −2 mol/m 3 was set close to upper estimate of this parameter, found in literature (Martinez-Cruz et al., 2015).
Due to high oxygen content, Kuivajärvi Lake is generally poor in methane (Miettinen et al., 2015).
To better understand the reasons for low surface methane concentrations, it is instructive to scrutinize 860 the budget of methane in the mixed layer (Fig. 14). Wee see that the methane fluxes nearly compensate each other, bubble fluxes being dominant in magnitude (see Table 3). Divergence of bubble flux is almost compensated by oxidation, whereas diffusion through thermocline is the smallest flux.
Thus, in the model, the epilimnetic and hypolimnetic pools of methane are almost "disconnected" due to minimal TKE in metalimnion (see Section 5.2). Moreover, the total methane influx from shal- In the numerical experiment "SS-" with LAKE model, where surface seiches (horizontal pressure gradient) were neglected, the seasonal pattern of methane concentration took the form presented at Fig. 9c. In this case, the basal methane content began to rise about 2 months earlier than it was ob-875 served (Fig. 9b) and calculated in the reference run (Fig. 9a) and reached maximal value of 598.5 mg/l vs. 351.5 mg/l in baseline experiment. This is caused by earlier O 2 depletion (not shown) due to negligible oxygen supply from above waters in conditions of very small TKE in hypolimnion (Fig.   6). Hence, we conclude that hypolimnetic turbulence is significant for gases accumulation and vertical distribution there, although it is likely to be of minor importance for mixed-layer concentrations 880 of O 2 , CO 2 and CH 4 , because of small gas transfer through metalimnion (see Section 5.2). are take into account. Standard k − turbulence closure is supplemented by parameterizations of internal seiches (Goudsmit, 2002), gravity waves (Mellor, 1989) and a new surface seiche formulation, following original concept by U.Svensson (Svensson, 1978).
The model is validated vs. extensive measurement data collected by University of Helsinki at Kuivajärvi Lake (Southern Finland) (Miettinen et al., 2015;Mammarella et al., 2015) during ice-905 free season of 2013 and including all meteorological variables above lake surface necessary to drive the model. In-water temperature, O 2 , CO 2 and CH 4 vertical profiles from the water column served to validate the model output.
The model was successful in capturing large-scale patterns of spatio-temporal variability of temperature and gases. In all the model parameterizations, only two constants relevant to CH 4 produc-910 tion and consumption were calibrated. The value of P 0 , regulating methane production in sediments, occurred to be very close to those obtained in our previous study for a thermokarst lake in North-Eastern Siberia (Stepanenko et al., 2011), corroborating the robustness of the model used. It is uncertainty in a number of other parameters, responsible for reactions involving O 2 and CO 2 , that is likely to contribute to model errors in hypolimnion and these of CO 2 in the surface layer.

915
As both carbon dioxide and methane typically accumulate below metalimnion in freshwater lakes (e.g. (Bastviken et al., 2008)), the vertical transport of these gases below mixed layer becomes an important factor for their evasion to the atmosphere. Our experiments together with stationary Richardson number analysis show that Mellor and Goudsmit extensions of k − model neither produce TKE in hypolimnion, nor generate turbulence in thermocline enough to sustain eddy diffusivity 920 above molecular constant. However, surface seiche parameterization allowed to produce turbulenceenhanced hypolimnion qualitatively consistent with empirical knowledge so far (Wüest et al., 2000;Wüest and Lorke, 2003). Reproducing considerable TKE in hypolimnion lead to much better correspondence of calculated CH 4 to observed one.
As there are strong doubts on complete suppression of turbulence even at Ri >> 1 (Saggio and so far, and should be addressed carefully in their future formulations. This will allow to get more rigorous regional and global estimates of greenhouse gases evasion to the atmosphere.      "GV+" -including gravity waves shear parameterization (Gill, 1982;Mellor, 1989), "IS+" -including internal seiche mixing parameterization (Goudsmit, 2002), "GV+SS-" -the same as "GV+" but with surface seiches switched off, "IS+SS-" -the same as "IS+" but with surface seiches switched off, "SS-" -the same as "base" but with