\documentclass[12pt,oneside]{book}
%aps,%
%12pt,%
%final,%
%notitlepage,%
%oneside,%
%onecolumn,%
%nobibnotes,%
%nofootinbib,% 
%superscriptaddress,%
%noshowpacs,%
%centertags]%
%{revtex4}

%\usepackage[T2A]{fontenc}
%\usepackage[koi8-r]{inputenc}
%\usepackage[russian]{babel}
\usepackage{amsmath}
\usepackage{authblk}
\usepackage{listings}
\usepackage{geometry}
\usepackage{array}
\usepackage{longtable}
\usepackage{gensymb}
\usepackage[font={Large}]{caption}
\usepackage{upgreek}
\usepackage[authoryear,round]{natbib}
\usepackage{graphicx}
\usepackage{hyperref}
%\usepackage[T1]{fontenc}

% Tables settings
\setlength\extrarowheight{4pt}

% Page geometry
 \geometry{
 a4paper,
 total={210mm,297mm},
 left=20mm,
 right=20mm,
 top=20mm,
 bottom=20mm,
 }

% Title
\title{The LAKE model. \\
Technical description}
\author{Victor Stepanenko}
\affil{Moscow State University}
\date{\today}


\begin{document}
 
\maketitle 

\chapter{The thermo- and hydrodynamics of the model}

\section{Basics of 1D equations in the water column}

Here we present a general 1D modelling framework used in the LAKE model. We start with the general transport Reynolds-averaged equation for the quantity $f$, that might be one of velocity components, temperature, turbulent kinetic energy (TKE) and other scalars:

\begin{equation}
\label{eq:fbasic}
\frac{\partial f}{\partial t} = - \frac{\partial u_i f}{\partial x_i} - \frac{\partial F_i}{\partial x_i} + R_f,
\end{equation}
assuming mass conservation equation for incompressible fluid:
\begin{equation}
\label{eq:contin}
\frac{\partial u_i}{\partial x_i}=0,
\end{equation}
where $u_i$ is the velocity component along $x_i$ Cartesian axis ($x_3=z$ being an axis pointing along gravity, $x_1=x,~x_2=y$), $F_i$ is the sum of non-advective fluxes of a property $f$ along $x_i$, and $R_f$ standing for the sum of sources and sinks of $f$. Now, introduce the averaging procedure as:

\begin{equation}
\overline{f}=\frac{\int_{A(z)}fdxdy}{A(z)},
\end{equation} 
with $A(z)$ denoting the horizontal cross-section of a lake. After applying this operator to \eqref{eq:fbasic} and making use of appropriate simplifications we get (for rigorous derivation we refer to Appendix \ref{app:aver}):

\subsection{Horizontal averaging}

\subsection{The normalized vertical coordinate}

\section{Thermodynamics of water column}

\section{Dynamics of water column}

\subsection{Governing equations}

\subsection{Parameterization of mean surface level gradient	}

In order to calculate the mean gradients of the surface level $\overline{\partial h_s / \partial x} , ~ \overline{\partial h_s / \partial y}$, we use the following parameterization, based on that originally proposed by U.Svensson \citep{Svensson1978}. 

First, for simplicity we assume the lake's body horizontal cross-section to be an ellipse at all levels. Then, we define South-North (S-N) and  West-East (W-E) vertical cross-sections of the lake's body, drawn through the ellipse center (given that the ellipses' centers at all levels fall onto the same vertical line, ref! to figure). Now, we denote the horizontal width of the S-N cross-section as $L_{S-N}(z)$, and define $L_{W-E}(z)$ analogously. The mass flux through the W-E section increases the lake surface level in the Northern half of a lake, and decreases it in the Southern part. The same holds for Western and Eastern parts. This means, that we can write

\begin{equation}
\frac{dh_N}{dt}A_0(t) = - \frac{dh_S}{dt}A_0(t) = 2\int_0^1 v L_{W-E} h d \xi,
\end{equation} 
\begin{equation}
\frac{dh_E}{dt}A_0(t) = - \frac{dh_W}{dt}A_0(t) = 2\int_0^1 u L_{S-N} h d \xi,
\end{equation}
where we introduced $h_S, h_N, h_W$ and $h_E$ as the average deviations of the lake surface level in Southern, Northern, Western and Eastern halfs of the lake surface, respectively, and $A_0(t)$ is the lake surface area. In the above equations we assumed that the deviations $h_*$ are small, so that the use of $A_0(t)$ instead of $A_0(h_*,t)$ is justified. 

Once $h_*$ are found, we parameterize the needed surface level gradients as:

\begin{equation}
g\overline{\frac{\partial h_s}{\partial x} } \approx \frac{g \pi^2}{4}  \frac{h_E - h_W}{L_{W-E,0}},
\end{equation}
\begin{equation}
g\overline{\frac{\partial h_s}{\partial y} } \approx \frac{g \pi^2}{4} \frac{h_N - h_S}{L_{S-N,0}}.
\end{equation}

It can be shown (ref!) that this parameterization of seiches exactly fits the analytical solution for specific case of 1D channel flow oscillations described by shallow water equations (Merian formula). To implement the above equations we need $L_{W-E}(z),L_{S-N}(z)$ for each particular lake to be simulated. These parameters are easily derived from $A(z)$, given the each lake horizontal cross-section is an ellipse, and its eccentricity is known (making the model cross-section shape to be close to the cross-section of a real lake).

\section{Turbulent mixing}

\section{Thermodynamics of ice}

\section{Thermodynamics of snow}

Ice cover acts in winter as a heat insulating layer and thereby controls the growth and melting of ice cover. Hence, the heat transport in a snow cover should be well reproduced. In order to acheive realistic results on snow temperature and depth, a number of processes are included in the model, snow gravitational compaction and liquid water transport, among others. The mathematical description of these processes closely follows the formulation from \citep{Volodina2000} (for liquid water equation derivation see Appendix \ref{app:snow}). The governing equations have a form:

\begin{equation}
\label{eq:snowtemp}
c_{sn} \rho_{sn} \frac{\partial T}{\partial t} = 
\frac{\partial}{\partial z} \lambda_{sn} \frac{\partial T}{\partial z} +
\rho_{sn}L_{fr}F_{fr},
\end{equation}
\begin{equation}
\label{eq:snowwat}
\frac{\partial W}{\partial t} = - (1-W) \frac{\rho_{w0}}{\rho_{sn}} \frac{\partial \gamma_v}{\partial z} - F_{fr} - \frac{W}{\rho_{sn}}C_{sn}.
\end{equation}
Here, the snow thermal conductivity is calculated by the empirical formula:

\begin{equation}
\lambda_{sn} = C_1 \left( \frac{\rho_{sn}}{\rho_{w0}} \right)^4 +
C_2 \frac{\rho_{sn}}{\rho_{w0}} + C_3,
\end{equation}
with $C_1=2.514,~C_2=0.796,~C_3=0.021$, $[C_{*}]=J/(m*s*K)$ \footnote{Hereafter, if not otherwise stated, * in the subscript means all possible values of the subscript in the current context}. The specific heat is determined by liquid water content:

\begin{equation}
c_{sn} = Wc_w + (1-W)c_i.
\end{equation}
The rate of phase transition, $F_{fr}$, is non-zero, when the temperature crosses the melting point for water, and is calculated according to heat balance in a given numerical cell. The gravitational flux of liquid water, m/s, is given by:

\begin{equation}
\gamma_v=h_c \left( \frac{W_v - W_{hc}}{p - W_{hc}} \right)^3,
\end{equation} 
where $W_v=W \rho_{sn}\rho_{w0}^{-1}$ is a volumetric ratio of liquid water. The snow compaction term, $C_{sn}$, enters the equation for $\rho_{sn}$:

\begin{equation}
\frac{\partial \rho_{sn}}{\partial t} = -\rho_{w0} \frac{\partial \gamma_v}{\partial z} + C_{sn}.
\end{equation}

Boundary condition at the top of snow is a heat balance equation, including heat and radiation fluxes. The botton boundary conditions, i.e. the condition at the snow-ice interface is continuity of heat flux and temperature. For liquid water content, the gravitational flux is zero at the bottom (boundary condition for the top is not needed, as \eqref{eq:snowwat} is a 1-st order equation).

\section{Thermodynamics of bottom sediments (ground)}

\chapter{Biochemistry in the water column and sediments}

\section{Governing equations for dissolved gases and organic carbon in a water column}

Thy dynamics of three gases is considered: methane ($CH_4$), oxygen ($O_2$) and carbon dioxide ($CO_2$). However, dissolved carbon dioxide is supposed to be always in carbonate equilibrium, so that it is concentration of dissolved inorganic carbon (DIC), $C_{DIC}=C_{CO_2}+C_{HCO_3^-}+C_{CO_3^{2-}}$, not of $CO_2$, that reflects the number of $CO_2$ molecules added to a solution (see Section \ref{sect:carbeq}). 

In addition, dissolved organic carbon (DOC), particulate organic carbon (both living, POCL, and dead, POCD) are calculated. 

The species listed above obey the following equation system:

\begin{align}
\label{eq:meth}
&\frac{\partial C_{CH_4}}{\partial t} = \text{Dif}_A(C_{CH_4}) + B_{CH_4} - O_{CH_4},\\
\label{eq:ox}
&\frac{\partial C_{O_2}}{\partial t} = \text{Dif}_A(C_{O_2}) + B_{O_2} + P_{O_2} - R_{O_2} - D_{O_2} - S_{O_2} - O_{O_2},\\
\label{eq:carbdi}
&\frac{\partial C_{DIC}}{\partial t} = \text{Dif}_A(C_{DIC}) + B_{CO_2} - P_{CO_2} + R_{CO_2} + D_{CO_2} + S_{CO_2} + O_{CO_2},\\
\label{eq:DOC}
&\frac{\partial \rho_{DOC}}{\partial t} = \text{Dif}(\rho_{DOC}) + E_{POCL} - D_{DOC},\\
\label{eq:POCL}
&\frac{\partial \rho_{POCL}}{\partial t} = \text{Dif}(\rho_{POCL}) + P_{POCL} - R_{POCL} - E_{POCL} - D_{h,POCL},\\
\label{eq:POCD}
&\frac{\partial \rho_{POCD}}{\partial t} = \text{Dif}(\rho_{POCD}) - \frac{w_g}{h}\frac{\partial \rho_{POCD}}{\partial \xi} - D_{POCD} + D_{h,POCL}.
\end{align}
where $\text{Dif}_A(\bullet)\equiv M(\xi,t) \frac{\partial \bullet}{\partial \xi} + \frac{1}{A h^2}\frac{\partial}{\partial \xi} \left( A k_s \frac{\partial \bullet}{\partial \xi} \right)$, $\text{Dif}(\bullet)\equiv M(\xi,t) \frac{\partial \bullet}{\partial \xi} + \frac{1}{h^2}\frac{\partial}{\partial \xi} \left( k_s \frac{\partial \bullet}{\partial \xi} \right)$, $M(\xi,t)=\left( \frac{\xi}{h} \frac{dh}{dt} - \frac{r-E}{h} \right)$ is a metric term arising from using the normalized vertical coordinate originating at moving water surface, $w_g$ is a sedimentation speed of POCD particles. Equations \eqref{eq:DOC}-\eqref{eq:POCD} do not contain $A$, i.e. they are not horizontally-averaged (see Appendix \ref{app:aver}) but ordinary 1D equations. This is caused by uncertainty how to assess marginal flux of these substances at the sediments-water interface. The r.h.s of these equations represent diffusion (assuming  $k_s$ to be the same eddy diffusivity for all species), sources and sinks due to the following processes:

\begin{itemize}
\item dissolution/exsolution of gases at the bubble-water interface ($B_{CH_4},~B_{O_2}$ and $B_{CO2}$);

\item photosynthesis ($P_{O_2},~P_{CO_2},~P_{POCL}$);

\item respiration ($R_{O_2},R_{CO_2},~R_{POCL}$);

\item biochemical oxygen demand in the water column ($D_{O_2},~D_{CO_2},~D_{DOC},~D_{POCD}$);

\item sedimentary oxygen demand ($S_{O_2},~S_{CO_2}$);

\item methane aerobic oxidation in the water column ($O_{CH_4},~O_{O_2},~O_{CO_2}$);

\item death of living species ($D_{h,POCL}$)
\end{itemize}

All variables in the above list are positive definite, excepting $B_{CH_4},B_{O_2}$ and $B_{CO2}$ that may be either positive or negative. All concentrations in \eqref{eq:meth}-\eqref{eq:carbdi} are expressed in $mol/m^3$ that allows for simple relations of sinks/sources in different equations based on stoichiometry of the respective reactions. In contrast, POC, DOCL and DOCD content in \eqref{eq:DOC}-\eqref{eq:POCD} is expressed as density, $kg/m^3$. Terms $B_{CO_2},~P_{CO_2},~R_{CO_2},~D_{CO_2},~S_{CO_2},~O_{CO_2}$ in \eqref{eq:carbdi} possess "$CO_2$" subscript because carbon atoms are supplied to or removed from DIC of a solution in a form of $CO_2$. 

We note that equations \eqref{eq:meth}-\eqref{eq:carbdi} and \eqref{eq:DOC}-\eqref{eq:POCD} may form a coupled system if sinks/sources at the r.h.s. are related, or be solved as independent subsets otherwise. The details will be given below.

In the following sections we first consider carbonate equilibrium, then continue with boundary conditions for \eqref{eq:meth}-\eqref{eq:carbdi} and finally describe the abovementioned sources/sinks in more detail. The formulations for photosynthesis, respiration, biochemical oxygen demand and sedimentary oxygen demand basically follow \citep{Stefan1994} and \citep{Hanson2004}.

\section{Carbonate equilibrium}\label{sect:carbeq}

Carbonate equilibrium means the equilibrium in the following reactions:

\begin{align}
&CO_2+H_2O \leftrightarrow H^+ + HCO_3^-,\\
&HCO_3^- \leftrightarrow H^+ +CO_3^{2-}. 
\end{align}
Involving kinetic constants of these reactions yields, that the DIC

\begin{equation}
\label{eq:carbeq}
C_{DIC}\equiv C_{CO_2}+C_{HCO_3^-}+C_{CO_3^{2-}} = C_{CO_2}\left[1+k_1 10^{pH}+k_1k_2 10^{2pH} \right].
\end{equation}
Here, the constants are given by Arrhenius equation:

\begin{equation}
k_i=k_{i0}\exp\left[-\frac{E_{act,i}}{R}\left(\frac{1}{T}-\frac{1}{T_0} \right) \right],~i=1,2,
\end{equation}
$R$ -- universal gas constant, $k_1=4.3* 10^{-7}~mol/l,~k_2=4.7* 10^{-11}~mol/l,~E_{act,1}=7.66* 10^3~J/mol,~E_{act,2}=1.49* 10^4~J/mol$. Thus, $C_{CO_2}$ is readily calculated given $C_{DIC}$ value, and vice versa. 

Carbon atoms are added or removed from carbonate equilibrium system in a form of $CO_2$ during respiration, photosynthesis and other processes, hence the change of $C_{DIC}$ equals to number of $CO_2$ consumed or produced. This explains equation \eqref{eq:carbdi}. For obtaining $CO_2$ flux across bubble surface or $CO_2$ diffusive flux to the atmosphere, $C_{CO_2}$ is calculated from \eqref{eq:carbeq}.

\section{Boundary conditions for dissolved gases in a water column}
The top boundary condition (at the lake-atmosphere interface) for any dissolved gas concentration for the case of open water has the form:

\begin{equation}
\left. \frac{k_s}{h}\frac{\partial C}{\partial \xi}\right	|_{\xi=0} = F_C,
\end{equation}

where $C$ is $C_{CH_4},C_{O_2}$ or $C_{CO_2}$, and $F_C$ is the diffusive flux of a gas into the atmosphere, positive upwards. This flux is calculated according to the widely used parameterization:

\begin{equation}
F_C = k_{ge}(C|_{\xi=0} - C_{ae}),
\end{equation}

with $C_{ae}$ being the concentration of the gas in water equilibrated with the atmospheric concentration and described by Henry law and $k_{ge},m/s$ denoting the gas exchange coefficient, the so-called "piston velocity". The latter is written as:

\begin{equation}
k_{ge} = k_{600}\sqrt{\frac{600}{Sc(T)}},
\end{equation}

with the Scmidt number $Sc(T)$ having individual values for different gases and being temperature-dependent (ref!). The $k_{600}$ coefficient has been a subject of numerous studies, and a number concepts have been proposed to quantify it \citep{Donelan2002}. 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 \citep{Cole1998}:

\begin{equation}
\label{eq:k_600}
k_{600} = C_{k_{600},1} + C_{k_{600},2} |\mathbf{u_{a,10}}|^{n_{k_{600}}}.
\end{equation}

Here, $\mathbf{u_{a,10}}$ stands for the wind speed at 10 m above the surface. The simple empirical eq. \eqref{eq:k_600} "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 on real lakes. Therefore, we also included surface renewal model \citep{MacIntyre2010,Heiskanen2014}, that in terms of $k_{600}$ states that:

\begin{equation}
\label{eq:k_600_SR}
k_{600}=\frac{C_{1,SR}(\epsilon|_{\xi=0} \nu_w)^{\frac{1}{4}}}{\sqrt{600}},
\end{equation}
where $\nu_w$ designates molecular viscosity of water. As TKE dissipation rate is available directly from $k-\epsilon$ closure, we do not use any special parameterization for $\epsilon|_{\xi=0}$ (e.g., \citep{MacIntyre2010}).

When a lake is covered by ice, $F_C=0$, neglecting contribution of diffusion through ice cracks. 

The boundary conditions for methane at the sediments-water body interface are:

\begin{equation}
-\left.\frac{k_s}{h}\frac{\partial C}{\partial \xi}\right|_{\xi=1}=-\left. k_{s,s}\frac{\partial C_s}{\partial z_s}\right|_{z_s=0},
\end{equation}
\begin{equation}
\label{eq:meth-lower_bc}
C|_{\xi=1}=\left. \left( \frac{C_s}{p} \right) \right|_{z_s=0}.
\end{equation}

Here, additional subscript "$s$" denotes characteristics of sediments. Porosity appears in \eqref{eq:meth-lower_bc} because $C_s$ is a bulk concentration in soil. These relations mean continuity of both diffusive flux and concentration at the water-sediments boundary. 

\begin{table}[h]
\caption{Constants for water-air gas exchange}
\begin{center}
\begin{tabular}{|l|c|c|c|}
\hline
Constant        & Variable in the code & Units & Value \\ \hline
$C_{k_{600},1}$ & constvel1            & $m/s$   & $5.75*10^{-6}$ \\ \hline
$C_{k_{600},2}$ & constvel2            & $(m/s)^{1-n_{k_{600}}}$      & $5.97*10^{-7}$ \\ \hline
$n_{k_{600}}  $ & wind10power          & $n/d$   & $1.7$          \\ \hline
$C_{1,SR}     $ & c1                   & $n/d$   & $0.5$          \\ \hline
\end{tabular}
\end{center}
\end{table}


\section{Dissolution/exsolution of gases at the bubble-water interface}

\subsection{Bubble model}
The bubble model presented here closely follows that described in \citep{mcginnis2006}.
Consider the evolution of a bubble, consisting of a gas mixture, rising from
a lake bottom. In this case the quantity of each $i$-th gas in the bubble $M_i$ (mol) 
is changing due to dissolution into oceanic water according to equation 

\begin{equation} \label{eq:mainbubble}
\frac{dM_i}{dt}=v_b\frac{dM_i}{dZ}=-4\pi r_b^2K_i(H_i(T)P_i-C_i), i=1,...,n_g, 
\end{equation}

where $r_b$ is the bubble radius (m), $H_i$ - the Henry constant  
dependent on temperature $T$ (K), $P_i$ the partial pressure of $i$-th gas (Pa), 
$C_i$ is the concentration of a gas dissolved in water
($mol/m^3$), $K_i$ is exchange coefficient, $v_b$ is bubble vertical velocity (m/s), $Z$ is the vertical coordinate originating at the bottom and pointing opposite to gravity (m), $n_g$ is the number of gases in a mixture. Five gases are considered simultaneously in a bubble:
methane, carbon dioxide, oxygen, nitrogen and argon. The temperature in the bubble is assumed to be equal to temperature of environment lake water at the depth of current bubble location. 
It means that the heat exchange between the rising bubble and water
is assumed to be intensive enough to dominate over the adiabatic cooling of the bubble. This is likely to be a rough aproximation, however, including adiabatic bubble temperature change would increase the complexity of the model. The temperature dependency of Henry constants for flat solute surface is taken from \citep{Sander1999}. The effect of gas-water interface curvature on equilibrium gas pressure is omitted because when using Thomson formula it turns to be negligible as long as typical bubble radiuses in lakes and ocean are considered ($\geq 1~mm$). Exchange coefficient $K_i$ is dependent on molecular diffusivity in water, bubble radius and its velocity according to empirical formulae from \citep{Zheng2002}. The bubble velocity is determined from equilibrium between buoyancy force and environment resistance defined by quadratic law for small radiuses ($r_b < 1.3 mm$) and taking into account bubble surface oscillations for larger sizes \citep{Jamialahmadi1994}.
\\ For each component of gas mixture one applies an ideal gas law because under the typical pressures at small depths (e.g. several dozens of meters) Van der Waals forces are small:

\begin{equation} \label{eq:gaslaw}
\frac{4}{3}P_i\pi r_b^3=M_iRT, i=1,...,n_g,
\end{equation}

where $R$ is the universal gas constant. The surface tension pressure is small for the bubbles with radiuses typical for lacustrine and marine environments. Then, when equating the gas mixture pressure $\sum_{i=1}^{n_g}P_i$ to hydrostatic pressure at a given depth $p_a+\rho_{w0}g(h-Z)$ ($\rho_{w0}$ being the mean density of water and $p_a$ the atmospheric pressure) and using \eqref{eq:gaslaw} one yields:

\begin{equation} \label{eq:r_b}
r_b=\left[\frac{3RT\sum_{i=1}^{n_g}M_i}{4\pi(p_a+\rho_{w0}g(h-Z))}\right]^{1/3}.
\end{equation}

For solution of $2n_g+1$ equations \eqref{eq:mainbubble}-\eqref{eq:r_b} the boundary conditions are needed. These are initial gases' quantities $M_{i,Z=0}=M_{i0}(t),i=1,...,n_g$, that are the quantities when the bubble crosses the lake bottom. In the model they are initialized as follows:

\begin{equation} \label{eq:initialbubble}
\begin{split}
\cr M_0=\frac{\frac{4}{3}\pi r_{b0}^3(p_a+\rho_{w0}gh)}{R\left.T\right|_{Z=0}},
\cr M_{i0}=\alpha_iM_0,i=1,...,n_g,
\end{split}
\end{equation}

where $M_0$ - is the total gas quantity in the bubble (mols). From \eqref{eq:initialbubble} it is seen that bubble initialization is given by inital bubble radius $r_{b0}$ and molar fractions of mixture components $\alpha_i$.
\\ The bubble model described above is numerically solved by Euler explicit scheme.

\subsection{Interaction between bubbles and dissolved gases}

To calculate $B_i,~B_1=B_{CH_4},~B_2=B_{CO_2},~B_3=B_{O_2}$ (see \eqref{eq:meth}-\eqref{eq:carbdi}) in our model we consider an idealized situation whenall bubbles rising from the bottom have the same initial radius $r_{b0}$ and identical gas composition. From \eqref{eq:mainbubble}-\eqref{eq:r_b} it immediately follows that their radius and composition are the same at any depth $\xi$. Since the vertical gas transport 
by bubbles changes with depth only due to exchange with water solution, one
writes

\begin{equation} \label{eq:biplus}
B_i=+\frac{1}{h}\frac{\partial F_{B,i}}{\partial \xi}.
\end{equation}

Here $F_{B,i}$ is the bubble gas flux ($mol/(m^2 s)$) pointed upwards (this leads to "+" sign in the r.h.s.). The definition of this flux is

\begin{equation}
F_{B,i}=M_i n_b v_b.
\end{equation}

We introduced bubble concentration ($m^{-3}$) here. All bubbles release at the bottom
and completely dissolve simultaneously at some depth or reach the surface. Moreover,
it is known \citep{Yamamoto2009,mcginnis2006} that bubbles with diameter
$\approx 1~mm$ are unstable and split up. Hence, in the model is it assumed that 
a bubble with $r_b \geq 0.5 cm$ splits into two. In the depth interval between two
subsequent bubble splits 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. One may rewrite \eqref{eq:biplus} as follows

\begin{equation}
B_i=\frac{1}{h}\frac{\partial (n_b v_b M_i)}{\partial \xi}=
\frac{F_{B0,i}}{h}\frac{\partial (Nm_i)}{\partial \xi},
\end{equation}

where we have introduced the relative bubble flux $Nm_i$ with $m_i=\frac{M_i}{M_{i0}}$,
$N=\frac{n_b v_b}{n_{b0} v_{b0}}$, and $F_{B0,i}$ is the bottom bubble flux (zero subscript indicates the bottom value of a variable). Evidently, $N(\xi)=2^k$, $k$ being the number of bubble divisions happened below the depth $\xi$. If the bottom bubble flux of one gas is known (in this model it is methane, $i=1$) then the bottom fluxes of other gases are determined by bottom bubble composition

\begin{equation}
F_{B0,i}=F_{B0,1}\frac{\alpha_i}{\alpha_1},i=2,...,n_g.
\end{equation}


\section{Photosynthesis}
The intensity of photosynthesis in terms of oxygen production is expressed as (ref!):

\begin{equation}
P_{O_2} = \dfrac{P_{max}L_{min}\rho_{Chl-a}}{H_{sec}\mu_{O_2}}. 
\end{equation}

The denominator here serves co convert the r.h.s. from $mg/(l*h)$ to $mol/(m^3 s)$. The $P_{max}$ value expresses limitation of oxygen production by temperature in a form:

\begin{equation}
P_{max} = C_{P}\theta_P^{(T-T_0)},
\end{equation}

so that $C_P$ is $P_{max}$ at the reference temperature $T=T_0$. The limitation of oxygen production by the available photosynthetically active radiation (PAR) is given by the following Haldane kinetics:

\begin{equation}
\label{eq:Lmin}
L_{min} = \frac{S_{PAR}(1 + 2 \sqrt{C_{Lmin,1}/C_{Lmin,2}})}
               {S_{PAR} + C_{Lmin,1} + S_{PAR}^2/C_{Lmin,2}} 
\end{equation}

The coefficients of this equation are specified as \citep{Stefan1994,Megard01011984}: 

\begin{equation}
C_{Lmin,1} = C_{PAR}\theta_{PAR}^{(T-T_0)},
\end{equation}
\begin{equation}
C_{Lmin,2} = \text{H}(T-T_{00})C_{Lmin,2,>T_{00}} + [1 - \text{H}(T-T_{00})]C_{Lmin,2,<T_{00}},
\end{equation}
with $\text{H}(\bullet)$ being a Heavyside function, and $T_{00}$ is an another reference temperature. It is seen from \eqref{eq:Lmin}, that $L_{min}\rightarrow0$ if $S_{PAR}\rightarrow 0$ and $S_{PAR}\rightarrow \infty$, i.e. PAR ihnibits photosynthesis at both low and high values of its intensity. To calculate $S_{PAR}$, it is assumed that the ratio of PAR intensity to total shortwave radiation intensity in the water column, $\alpha_{PAR}$, is constant with depth, thereby

\begin{equation}
S_{PAR}=\alpha_{PAR}H_{sec}T_{J\rightarrow Eins}S.
\end{equation}  

The coefficient transfroming from $J$ to $Einstein$ (Einstein is an energy of Avogadro number of photons), $T_{J\rightarrow Eins}$, is calculated in approximation of uniform distribution of energy in PAR region, yielding

\begin{equation}
T_{J\rightarrow Eins} = \frac{\lambda_{PAR}}{N_A h_P c},
\end{equation}

with $N_A, h_P, c$ denoting Avogadro number, Planck constant and the light speed in vacuum, respectively, all in SI units.

Finally, from the photosynthesis reaction:

\begin{equation}
\label{eq:photosynt}
6CO_2 + 12H_2O + photons \rightarrow C_6H_{12}O_6 + 6O_2 + 6H_2O,
\end{equation}
or, in a shortened form:
 
\begin{equation}
CO_2 + 2H_2O + photons \rightarrow CH_2O + O_2 + H_2O,
\end{equation}

we see that the carbon dioxide consumption equals oxygen production, i.e. $P_{CO_2} = P_{O_2}$.

\begin{table}[h]
\caption{Constants in photosynthesis model}
\begin{center}
\begin{tabular}{|l|c|c|c|}
\hline
Constant        & Variable in the code  & Units    & Value \\ \hline
$C_P$           & c1\_pmax              & $h^{-1}$ & 9.6   \\ \hline
$\theta_P$      & c2\_pmax              & n/d      & 1.036 \\ \hline
$T_0$           & T\_0                  & $\celsius$ & 20  \\ \hline
$T_{00}$        & T\_00                 & $\celsius$ & 10  \\ \hline
$\mu_{O_2}$     & molmass\_o2           & $g/mol$  & 32    \\ \hline
$H_{sec}$       & hour\_sec             & $s$      & 3600 \\ \hline
$C_{PAR}$       & k1\_c1                & $Einstein/(m^2*h)$ & 0.687 \\ \hline
$\theta_{PAR}$  & k1\_c2                & n/d      & 1.086 \\ \hline
$C_{Lmin,2,>T_{00}}$ & k2\_c2           & $Einstein/(m^2*h)$ & 15. \\ \hline
$C_{Lmin,2,<T_{00}}$ & k2\_c1           & $Einstein/(m^2*h)$ & 5. \\ \hline
$\lambda_{PAR}$ & lambda\_PAR0          & m        & $5.5*10^{-7}$ (550 nm) \\ \hline
$\alpha_{PAR}$  & short2PAR             & n/d      & 0.48 \\ \hline
\end{tabular}
\end{center}
\end{table}

Now, we convert photosynthesis rate to appropriate units of $P_{POCL}$. This term is a mass gain of $CH_2O$ in living species, hence, using the above stoichiometry, we have:

\begin{equation}
P_{POCL} = P_{CO_2}\mu_{CH_2O}.
\end{equation} 

\section{Respiration}

\subsection{Formulation by Stefan and Fang}

Respiration is a process opposite to photosynthesis. In the model, the oxygen consumption due to respiration is given by:

\begin{equation}
R_{O_2} = \frac{k_r\theta_r^{T-T_0}\rho_{Chl-a}}{Y_{CHO2}D_{sec}\mu_{O_2}},
\end{equation}

where $Y_{CHO2}$ is a mass ratio of chlorophyll-a to oxygen, utilized in respiration. Analogously to the case of photosynthesis, $R_{CO_2} = R_{O_2}$.

\begin{table}[h]
\caption{Constants in respiration model}
\begin{center}
\begin{tabular}{|l|c|c|c|}
\hline
%\title{List of constants for water-air gas exchange}\\ \hline
Constant        & Variable in the code  & Units    & Value     \\ \hline
$Y_{CHO2}$      & YCHO2                 & $mg/mg$  & $8*10^{-3}$ \\ \hline
$k_r$           & k\_r                  & $day^{-1}$ & $1*10^{-1}$ \\ \hline
$\theta_r$      & theta\_r              & n/d       & 1.045 \\ \hline
$D_{sec}$       & day\_sec               & s         & 86400 \\ \hline
$\mu_{O_2}$     & molmass\_o2           & $g/mol$  & 32    \\ \hline
\end{tabular}
\end{center}
\end{table}

\subsection{Formulation by Hanson et al.}

P.Hanson et al. \citep{Hanson2004} assume, that respiration is performed by "living particles", POCL, only in epilimnion, and may be scaled by gross primary production (i.e., photosynthesis rate),  $R_{POCL}=\alpha_{POCL}P_{POCL},~\alpha_{POCL}=0.8$.  
In contrast, we assume this process to happen unless not enough oxygen \emph{in situ} is available. The value of $R_{POCL}$ is a mass loss of organic matter, that can be generally represented as $CH_2O$, and from \eqref{eq:photosynt} one sees that the respective quantity of $CO_2$ produced is $R_{POCL}/\mu_{CH_2O},~mol/(m^3*s)$. Thus,

\begin{equation}
R_{O_2}=R_{CO_2}=\frac{\alpha_{POCL}P_{POCL}}{\mu_{CH_2O}}.
\end{equation}

\section{Biochemical oxygen demand (BOD)}

\subsection{Formulation by Stefan and Fang}

Biochemical oxygen demand is an oxygen consumption due to organic matter decomposition in the water column. It is expressed as \citep{Stefan1994}:

\begin{equation}
B_{O_2} = \frac{k_b \theta_b ^{T-T_0}\rho_{Chl-a}C_{st}}{D_{sec}\mu_{O_2}}. 
\end{equation}

Here, temperature dependence constant is given by:

\begin{equation}
\theta_b = [\theta_{b1}\text{H}(T-T_0) + \theta_{b2}[1-\text{H}(T-T_0)]]\text{H}(T-T_{md}).
\end{equation}

This means switching between constants $\theta_{b1}$ and $\theta_{b2}$ at the reference temperature $T_0$ and "switching off" BOD below the temperature of maximum density, $T_{md}$, implying ice-cover conditions. The constant $C_{st}$ reflects the stoichiometry of respective reactions. According to \citep{Stefan1994}, it reads:

\begin{equation}
C_{st} = M_{O_2/C}M_{C/Chl},
\end{equation}

where $M_{O_2/C}$ is the mass ratio of oxygen to carbon in the aerobic organic decay reaction:

\begin{equation}
\label{eq:BOD}
C + O_2 = CO_2,
\end{equation}

and $M_{C/Chl}$ is a mass ratio of carbon to chorophyll-a in the organic matter. The equation \eqref{eq:BOD} shows that $D_{CO_2} = D_{O_2}$.

\begin{table}[h]
\caption{Constants in biochemical oxygen demand model}
\begin{center}
\begin{tabular}{|l|c|c|c|}
\hline

Constant        & Variable in the code  & Units    & Value     \\ \hline
$k_b$           & k\_b                  & $day^{-1}$ & 0.1 \\ \hline
$\theta_{b1}$     & theta\_b1             & n/d        & 1.047 \\ \hline
$\theta_{b2}$     & theta\_b2             & n/d        & 1.13 \\ \hline
$T_{md}$          & T\_md                 & \celsius   & 4. \\ \hline
$M_{O_2/C}$       & mO2C\_dec              & mg/mg      & 8/3 \\ \hline
$M_{C/Chl}$       & mCChla\_dec            & mg/mg      & 30. \\ \hline
\end{tabular}
\end{center}
\end{table}

\subsection{Formulation by Hanson et al.}

Here, for biochemical oxygen demand, we adopt respiration terms, related to dead organic particles (POCD) and dissolved organic carbon (DOC) from \citep{Hanson2004}, namely, $D_{POCD}$ and $D_{DOC}$ in our notation. P.Hanson et al. suggest that $D_{POCD}=\rho_{POCD}/\tau_{POCD},
~D_{DOC}=\rho_{DOC}/\tau_{DOC}$ with time scales $\tau_{POCD}=20D_{sec},~\tau_{DOC}=
200D_{sec}$. Again, this is a mass loss by POCD and DOC, and the conversion to molar concentration change rate leads to:

\begin{equation}
D_{O_2}=D_{CO_2}=\frac{1}{\mu_{CH_2O}}\left(\frac{\rho_{POCD}}{\tau_{POCD}} + \frac{\rho_{DOC}}{\tau_{DOC}} \right).
\end{equation}


\section{Sedimentary oxygen demand (SOD)}
Two sedimentary oxygen demand formulations are implemented in the model. In both of them, the sedimentary oxygen demand appears as a sink in \eqref{eq:ox} and in essence is the contribution of the vertical flux at the lake's bottom to the horizontally averaged oxygen concentration:

\begin{equation}
S_{O_2} = - \frac{F_{SOD}}{Ah}\frac{\partial A}{\partial \xi}.
\end{equation}

In the first SOD formulation which is that used in \citep{Stefan1994}, the bottom oxygen flux due to organics decomposition is a function of temperature only:

\begin{equation}
  F_{SOD} = \frac{1}{D_{sec}\mu_{O_2}10^3} S_{b20}\theta_s^{T-T_0}.
\end{equation}
Here, again switching between two different values of temperature dependence parameter, $\theta_s$,  is used:

\begin{equation}
  \theta_s = \theta_{s1}\text{H}(T-T_0) + \theta_{s2}[1-\text{H}(T-T_0)].
\end{equation}

This parameterization does not include dependence of SOD on oxygen concentration, hence the $O_2$ uptake in sediments is not zero even when $C_{O_2}=0$. From this point the more physically grounded SOD parameterization is the model provided by Robert Walker and William Snodgrass  \citep{Walker1986}. Basing on the argument that SOD is contributed both by diffusion (governed by Fickian law) and biochemical consumption (described by Michaelis-Menthen kinetics), they derive:

\begin{equation}
  F_{SOD} = \mu_{\beta}\frac{C_{O_2}}{k_{O_2,SOD}+C_{O_2}} + k_c C_{O_2},
\end{equation}
where $\mu_{\beta}$ is proportional to organics oxidation potential rate in sediments, and $k_c$ is the mass transfer coefficient. Both are thought to be exponentially dependent on temperature:

\begin{equation}
  \mu_{\beta} = \mu_{\beta,0}\theta_{\mu_{\beta}}^{T-T_{\mu_{\beta}}},
  k_c = k_{c,0}\theta_{k_c}^{T-T_{k_c}}.
\end{equation}

The stoichiometry of SOD is assumed to be close to that of BOD \eqref{eq:BOD}, therefore, $S_{CO_2} = S_{O_2}$. Additionally, the flux of $O_2$ due to SOD at the lake bottom, $F_{SOD}$, is used as the bottom (lake deepest point) boundary condition for the oxygen equation \eqref{eq:ox}.

\begin{table}[h]
\caption{Constants in sedimentary oxygen demand model}
\begin{center}
\begin{tabular}{|l|c|c|c|}
\hline

Constant        & Variable in the code  & Units    & Value     \\ \hline
$S_{b20}$       & Sb20                  & $mg/(m^2day)$ & $10^3$ \\ \hline
$\theta_{s1}$     & theta\_s1             & n/d        & 1.065 \\ \hline
$\theta_{s2}$     & theta\_s2             & n/d        & 1.13 \\ \hline
$\theta_{\mu_{\beta}}$  & thetamu\_SOD    & n/d        & 1.085 \\ \hline
$\theta_{k_c}$     & thetaC\_SOD          & n/d        & 1.103 \\ \hline
$T_{\mu_{\beta}}$  & -             & K        & 25 \\ \hline
$T_{k_c}$          & -             & K        & 20 \\ \hline
$\mu_{\beta,0}$    & mubeta0       & $mol/(m^2*s)$ & $0.5/(\mu_{O_2}D_{sec}),~[\mu_{O_2}]=g/mol$ \\ \hline
$k_{c,0}$          & kc0           & m/s        & 0.045/$D_{sec}$ \\ \hline
\end{tabular}
\end{center}
\end{table}

\section{Exudates and death rate of POCL}

Hanson et al. suggest exudation to be scaled with photosynthesis rate, $E_{POCL}=\beta_{POCL}P_{POCL},~\beta_{POCL}=0.03$
 and the death rate to be defined as $D_{h,POCL}=\frac{\rho_{POCL}}{\tau_{Dh}}$, where time scale $\tau_{Dh}$ ranges from $1.1D_{sec}$ in hypolimnion to $33D_{sec}$ in epilimnion.
 
 
\section{Sedimentation of organic particles}

In the current model version we use the Stokes sedimentation velocity below the mixed layer

\begin{equation}
w_s=\frac{4}{3A}\frac{\Delta g d^2}{\nu_m},
\end{equation}
and the high-Reynolds-number limit of this variable

\begin{equation}
w_s=\sqrt{\frac{4}{3B}\Delta g d}
\end{equation}
in the mixed layer. Here, $\Delta=\rho_{p}/\rho_{w0}-1$, $\rho_p$ is a particle's density, and $d$ -- its diameter,  the typical values for constants may be chosen as $A=30.0$, and $B=1.25$  \citep{Song2008}, and the density of organic particles as $1.25~g/cm^3$ \citep{Avnimelech2001}. 


\section{Chlorophyll-a dynamics}

The chlorophyll-a dynamics in the model follows a simple scheme suggested in \citep{Stefan1994}. In it, chlorophyll-a density is calculated as:

\begin{equation}
\rho_{Chl-a} = \rho_{Chl-a,0}\text{H}(H_a - z),
\end{equation}

where the active layer, $H_a$, is a maximum between mixed-layer depth, $H_{ML}$, and the photic zone depth, $H_{PZ}$. The mixed-layer depth is defined as the depth of maximum Brunt-V\"{a}is\"{a}l\"{a} frequency, and the photic zone depth is estimated as the depth at which the PAR irradiance is 10\% of its surface value. The mean chlorophyll-a concentration, $\rho_{Chl-a,0}$, is assigned according to a trophic status of the lake: $2*10^{-3}~mg/l$ for oligotrophic lakes, $6*10^{-3}~mg/l$ for mesotrophic lakes and $15*10^{-3}~mg/l$ for eutrophic lakes. In turn, the trophic status is formally defined from the water turbidity. The Secchi disk values of 2 m and 3.5 m are used to distinguish between eutrophic and mesotrophic, mesotrophic and oligotrophic states, respectively. These thresholds are expressed in the model through light extinction coefficient values, $\alpha$, using Poole and Atkins formula \citep{Poole2009}:

\begin{equation}
\alpha = \frac{k_{PA}}{z_{SD}}, 
\end{equation}

where $z_{SD}$ is the Secchi disk depth and $k_{PA}=1.7$. The above chlorophyll-a scheme is identical to that of \citep{Stefan1994}, excepting for it does not take into account the annual cycle of $\rho_{Chl-a,0}$.



\section{Aerobic methane oxidation}
The aerobic oxidation of methane follows the equation:

\begin{equation}
CH_4 + 2O_2 = CO_2 + 2H_2O
\end{equation}

that means $O_{O_2} = 2O_{CO_2} = 2O_{CH_4}$. To calculate $O_{CH_4}$, three options are available in the model. The simplest one is the first-order kinetics:

\begin{equation}
O_{CH_4} = k_{ox,CH_4}C_{CH_4},
\end{equation}

where $k_{ox,CH_4}$ is a constant (according to \citep{Striegl1998}, is chosen $0.38D_{sec}^{-1}s^{-1}$). 

More sophisticated, but containing 4 constants, Michaelis-Menten kinetics, has the general form:

\begin{equation}
\label{eq:M-M_methox}
O_{CH_4} = V_{max} exp \left[ \frac{\Updelta E_{ox,CH_4}}{R} 
           \left( \frac{1}{T} - \frac{1}{T_0} \right) \right] 
           \frac{ C_{CH_4}}{K_{hs,CH_4} + C_{CH_4}} 
           \frac{ C_{O_2}}{K_{hs,O_2} + C_{O_2}},
\end{equation}

where we introduce the activation energy of methane oxidation reaction, $\Updelta E_{ox,CH_4}$, reaction potential $V_{max}$ at the reference temperature, $T_0$, half-saturation constants $K_{hs,CH_4}, K_{hs,O_2}$. However, a number of studies \citep{Utsumi1998,Liikanen2002,Striegl1998} have shown a weak dependence of reaction potential on temperature. Hence, the second option in the model is the Michaelis-Menten kinetics \eqref{eq:M-M_methox} without exponential term.

In a case of oxygen depletion aerobic methane oxidation should approach zero in the water column, so that $O_2$ should be kept in \eqref{eq:M-M_methox}. Unfortunately, we have found a quite limited data on $K_{hs,O_2}$ in the literature. However, from the experience on methane modeling in wetlands \citep{Watson1997} and measurements in the lake sediments \citep{Lidstrom1984} $K_{hs,O_2}\approx K_{hs,CH_4}$, and we have used this approximation in our model.

In order to avoid dependency of methane oxidation on oxygen concentration the third option considers simplified Michaelis-Menten equation, formally coming from \eqref{eq:M-M_methox} for highly oxygenated waters:

\begin{equation}
O_{CH_4} = V_{max} \frac{ C_{CH_4}}{K_{hs,CH_4} + C_{CH_4}}.
\end{equation}

The Michaelis-Menten constants are given in the table \ref{table:M-M-methoxid}.

\begin{table}[h]
\caption{Michaelis-Menten kinetics constants for methane oxidation in the water column}
\label{table:M-M-methoxid}
\begin{center}
\begin{tabular}{|lp{4cm}|cp{3cm}|cp{3cm}|cp{3cm}|}
\hline

Reference        & $V_{max}$  & $K_{hs,CH_4}$ & Remarks  \\ \hline
\citep{Lidstrom1984}   & $3.8*10^{-2}H_{sec}^{-1}~mol/(m^3s)$  & $9.5\pm1.2~mol/m^3$ & Measured in top 1 cm of bottom sediments \\ \hline
\citep{Liikanen2002}   & $3.6*10^{-2}D_{sec}^{-1}~mol/(m^3s)$   & $5.5*10^{-3}~mol/m^3$  & shallow profundal water, 4 m  \\ \hline
\citep{Liikanen2002}   & $1.4*10^{-1}D_{sec}^{-1}~mol/(m^3s)$   &   $4.4*10^{-2}~mol/m^3$  & deep profundal water, 9 m \\ \hline
\end{tabular}
\end{center}
\end{table}

\section{Methane dynamics in the bottom sediments}

\subsection{Governing equation}

The dynamics of bulk methane concentration (the number of dissolved molecules in moles per unit volume of sediment / soil) is governed by three processes: production, ebullition and diffusion:

\begin{equation}
\label{eq:methsoil}
\frac{\partial C_{CH_4}}{\partial t} = \frac{\partial}{\partial z_s} \left(k_{CH_4} \frac{\partial C_{CH_4}}{\partial z_s}  \right) + P_{soil,CH_4} - E_{soil,CH_4}- O_{soil,CH_4},
\end{equation}

where $z_s$ is a downward coordinate originating at the sediments' surface, $P_{soil,CH_4}$ and $E_{soil,CH_4}$ are the methane production and ebullition rates, respectively, and $O_{soil,CH_4}$ is the aerobic methane oxidation, all positive. Aerobic methane oxidation term is non-zero only in the top thin computational layer; in the rest layers oxygen is assumed to be depleted, that is the well-known observational fact \citep{Huttunen2006}. Note, that the methane concentration $C_{CH_4}$ accounts for dissolved gas only, so that ebullition rate is a sink term for it. The soil diffusion coefficient is given by:

\begin{equation}
k_{CH_4} = C_{tort} k_{CH_4,w},
\end{equation}

and molecular diffusivity for methane dissolved in water, $k_{CH_4,w}$, depends on temperature according to quadratic function. Here, $C_{tort}$ is a tortuosity coefficient (set to 0.66), accounting for tortuosity of diffusion paths through soil skeleton.

The eq. \eqref{eq:methsoil} should also contain the plant-mediated transport of methane for vegetated lakes. In future we intend to include this mechanism in the model.

\subsection{Production}

Methane production is comprised of two sources \citep{Stepanenko2011}. First is methane production from decomposition of "young" organics, i.e. the current NPP of lake ecosystem, depositing at the lake bottom. The organic material reaching the bottom after coast abrasion also contributes to this source, located at the top of sediments. The second source of methane is to be switched on in the model for thermokarst lakes for which it is known, that the "old" organics, sequestered in the premafrost below the talik, is subjected to biotic degradation when the temperature exceedes the metling point in the process of talik deepening \citep{Walter15072007}. Thus, this second source has its maximum in the vicinity of talik bottom. 

Formulations for both production sources base at the assumption that the production is proportional to quantity and quality of organics, exponentially dependent on temperature and does not happen under temperatures less than melting point:

\begin{equation}
\label{eq:CH4prodgen}
P_{soil,CH_4,i} \propto \rho_{org,i}\text{H}(T-T_{mp})q_{10}^{T/10}(1+\alpha_{O_2,inhib}C_{O_2})^{-1}
\left[1-\text{H}(C_{SO_4}-C_{SO_4,crit})\right],~~~i=old,new
\end{equation}
where $T_{mp}$ standing for the melting point temperature, $q_{10}$ is times that the methane production increases when the temperature rises by $10 \celsius$, and $\alpha_{O_2,inhib}$ denoting the constant for methanogenesis inhibition by oxygen. Oxygen concentration, $C_{O_2}$, decays very fast within a few millimeters below the bottom \citep{Huttunen2006}. The last multiplier in \eqref{eq:CH4prodgen} represents suppression of methanogenic \emph{Archaea} by sulfate-reducing bacteria when sulfate concentration exceeds critical value, $C_{SO_4,crit}$ \citep{Lovley1986}. The constant of proportionality in \eqref{eq:CH4prodgen} should reflect the quality of the organic substrate regarding methane production. The density ($kg/m^3$) of organics available for anaerobic decomposition resulting in methane production, $\rho_{org,i}$, may be estimated from either $1^{st}$-order or Michaelis-Menten kinetics. For young organics we adopt an assumption $\rho_{org,new} \propto \exp(-\alpha_{new}z_s)$ \citep{Walter2000}, that formally can be derived from $1^{st}$-order kinetics for $\rho_{org,new}$ decay and constant rate of depositing new sediments at the lake bottom. For the old organics content, an expression is developed in \citep{Stepanenko2011} based on Michaelis-Menten decomposition rate for $\rho_{org,old}$:

\begin{equation}
\label{eq:oldorgdens}
\rho_{org,old}=\rho_{org,old,0} \left[ 2 + \lambda_{\rho} -
\sqrt{(1 + \lambda_{\rho})^2 + \gamma_{\rho} (h_t^2 - z_s^2)} \right],
\end{equation}
 
where $\lambda_{\rho} = {\rho}_{org,old,0}K_{org,old}^{-1},~\gamma_{\rho} = 2 V K_{org,old}^{-1} C_t^{-2}$ are the derived constants, and the original constants are given in the Table \ref{table:Methgen_sed}, $h_t$ - talik depth. Thus, basing on \eqref{eq:oldorgdens}, we get:

%\begin{equation}
\begin{align}
P_{soil,CH_4,i} = P_{i,0} \rho_i^* \text{H}(T-T_{mp})q_{10}^{T/10}(1+\alpha_{O_2,inhib}C_{O_2})^{-1}
\end{align}
\[\rho_i^* = \left\{
  \begin{array}{lr}
    \exp(-\alpha_{new}z_s) & : i=new\\
    \left[ 2 + \lambda_{\rho} -
\sqrt{(1 + \lambda_{\rho})^2 + \gamma_{\rho} (h_t^2 - z_s^2)} \right] & : i=old
  \end{array}
\right.
\]
%\end{equation}
Here, $P_{new,0},P_{old,0}$ are new constants to be calibrated. 
 
\begin{table}[h]
\caption{Constants for methane production in the lake sediments}
\label{table:Methgen_sed}
\begin{center}
\begin{tabular}{|p{4cm}|p{3cm}|p{3cm}|p{3cm}|}
\hline

Constant        & Variable in the code  & Units    & Value (reference)    \\ \hline
${\rho}_{org,old,0}$ (the density of old organics below the talik)    
& C0\_oldorg & $kg/m^3$ & 18 \\ \hline
$K_{org,old}$ (the half-saturation constant of Michaelis-Menten decomposition of old organics in talik)    
& k\_oldorg             & $kg/m^3$        & $3*10^{-1}$ \\ \hline
$V$ (the maximum decomposition rate of old organics in talik)     
& V\_oldorg             & $kg/(m^3*yr)$        & $2*10^{-3}$ \\ \hline
$C_t$ (a constant in a talik deepening law, $h_t=C_t \sqrt{t}$, $t$ - time)     
& $ms^{-1/2}$           & C\_talik\_age        & 0.5 \citep{West2008} \\ \hline
$q_{10}$     
& q100         & n/d     & 2.3 \citep{Liikanen2002a} \\ \hline
$\alpha_{new}$     
& lambda\_new\_org       &  $m^{-1}$       & 3 \\ \hline

\end{tabular}
\end{center}
\end{table} 
 
\subsection{Ebullition}

The formation of bubbles (ebullition) happens in the model only in the sediments layer, as we havn't found evidences in literature that this process occurs at significant rates in the water column. As the bubble formation takes a certain amount of methane from the dissolved state, it is represented by a sink term in \eqref{eq:methsoil}. The rate of this sink is determined by a methane concentration excess over a critical value, $C_{CH_4,crit}$:

\begin{equation}
E_{soil,CH_4} = r_{ebul}\text{H} \left (C_{CH_4} - C_{CH_4,crit} \right )
\left(1-\text{H}(w_i) \right) \left(C_{CH_4} - C_{CH_4,crit} \right),
\end{equation} 
where we also postulated that ebullition happens in non-frozen soils only, and $r_{ebul}=1.0~H_{sec}^{-1}$ \citep{Walter2000}. Now, the critical concentration of dissolved methane in soil, $C_{CH_4,crit}$, is assumed to be a fraction, $r_{crc}=0.4$, of saturation concentration (\citep{Wania2007} and references therein). In turn, the saturation concentration is found from two conditions: (i) the sum of partial pressures of all gases comprising a bubble equilibrates the external hydrostatic load, $P=g(\rho_{w0}h+\rho_s z_s)$, and (ii) partial pressure of each gas is in equilibrium with the dissolved gas according to Henry's law (with Henry constants $H_*(T)$, * standing for the gas name). We assume the presence of two dissolved gases in sediments, entering the bubble composition: nitrogen and methane. Thus, the critical bulk methane concentration is expressed as:

\begin{equation}
\label{eq:ch4crit}
C_{CH_4,crit} = r_{crc}p \left( PH_{CH_4}(T) - \frac{H_{CH_4}(T)}{H_{N_2}(T) } C_{N_2} \right),
\end{equation}
with $p$ standing for sediments' porosity. Following \citep{Bazhin2001} we assume that nitrogen rapidly decays to zero in the top thin sediments layer, and therefore contributes to \eqref{eq:ch4crit} only in the top numerical layer of the in-sediments' grid.

The mechanical interaction of bubbles with the soil skeleton \citep{Scandella2011} is omitted in the current version of the model, and it is assumed that all bubbles generated at different depths in the bottom sediments instantly reach the lake's bottom. Therefore, the methane bubble flux at the lake bottom is:

\begin{equation}
F_{CH_4,b}(\xi=1) = \int_0^{d} E_{soil,CH_4} dz_s,
\end{equation}
where $d$ is the thickness of ground domain in the model. 

\subsection{Aerobic methane oxidation in bottom sediments}

\chapter{List of symbols}

\begin{center}
\begin{tabular}{ p{3cm} p{5cm} p{3cm} p{3cm} }
\hline
Symbol     & Name   & Variable in the code & Units   \\ \hline
\textbf{General} \\ \hline
$h$        & Lake depth     & h1, h2       & m       \\ \hline
$l$        & Ice  thickness & l1, l2       & m       \\ \hline
$h_{sn}$   & Snow thickness & hs1          & m       \\ \hline
$\xi=\frac{z}{h}$      & Normalized vertical coordinate & dzeta  & n/d       \\ \hline
\textbf{Snow cover}  \\ \hline
$W$        & Liquid water content in snow  & - & kg/kg       \\ \hline
\textbf{Biochemistry} \\ \hline
$\rho_{Chl-a}$  & Chlorophyll-a density in a water column  & Chl\_a  & $mg/l$  \\ \hline   
$S_{PAR}$  & The intensity of photosynthetically active radiation  & PAR  & $Einstein/(m^2*h)$  \\
\hline
\end{tabular}
\end{center}

\bibliographystyle{plainnat}
\bibliography{lake_methane&oxygen,wetland_methane,Stepanenko,misc}

\appendix

\chapter{Equation for horizontally averaged quantity in a lake}
\label{app:aver}

Consider equation \eqref{eq:fbasic} and an auxiliary operator:

\begin{equation}
\label{eq:auxav}
\tilde{f}=\int_{A(z)}fdxdy.
\end{equation}
The cross-section of a lake with notations used in this derivation is given at Fig. \ref{fig:A}.

\begin{figure}[h!]
\label{fig:A}
\centering
\includegraphics[scale=1]{Pic1.png}
\caption{A lake horizontal cross-section}
\end{figure}
The integration operator \eqref{eq:auxav} possesses the following property:

\begin{equation}
\frac{\partial \tilde{f}}{\partial z} = \widetilde{\frac{\partial f}{\partial z}} + B_f,
\end{equation}
\begin{equation}
B_f=\int_{x_1(z)}^{x_2(z)}\left[ \frac{\partial y_2}{\partial z}f(x,y_2,z) - \frac{\partial y_1}{\partial z}f(x,y_1,z) \right]dx,
\end{equation}
stemming from the Leibnitz integral rule. Now let's apply operator $\widetilde{(...)}$ to \eqref{eq:fbasic}, then insert $\overline{f}=A\tilde{f}$, and it will lead us to

\begin{equation}
\label{eq:fbasicav}
\frac{\partial A\overline{f}}{\partial t}=-\int_{\Gamma_{A(z)}}f(\mathbf{u_h}\cdot\mathbf{n})dl - 
\int_{\Gamma_{A(z)}}(\mathbf{F_h}\cdot\mathbf{n})dl - \frac{\partial A\overline{wf}}{\partial z} - 
\frac{\partial A\overline{F_z}}{\partial z} + B_{wf} + B_{F_z} + A\overline{R_f},
\end{equation}
where we introduced $\mathbf{u_h}=\{u_1,u_2\}$ and $\mathbf{F_h}=\{F_1,F_2\}$, and $\Gamma_{A(z)}$ is a boundary of $A$ at depth $z$. The first term to the right hand side of \eqref{eq:fbasicav} is a horizontal advection of property $f$ through boundaries of a water basin, i.e. the inflow from inlets, outflow by outlets and groundwater discharge. The second term represents non-advective horizontal fluxes at the lake margins, whereas $B_{*}$ quantifies the effect of vertical fluxes at the lake bottom of depth $z$; we refer to them as to \emph{marginal fluxes}. Equation \eqref{eq:fbasicav} is the most general equation, that is, however, difficult to implement without further simplifications. First, assume that the lake bottom is quasi-horizontal, and in this case the rigid boundary condition for velocity brings $w \approx 0,~B_{wf} \approx 0$. Then, we suppose that $\mathbf{F}=\{F_1,F_2,F_3\}$ is normal to the bottom boundary, and it is true for diffusive transport, because it is proportional to a gradient of $f$, and this gradient is usually oriented almost perpendicular to the bottom surface. Therefore, $F_1 \approx 0,~F_2 \approx 0$, vanishing the second term to the right hand side of \eqref{eq:fbasicav}. We also can decompose the vertical advection as $\overline{wf}=\overline{w}\overline{f}+\overline{w'f'},~w'=w-\overline{w},~f'=f-\overline{f}$. After these modifications, \eqref{eq:fbasicav} devolves to:

\begin{equation}
\label{eq:fbasicav2}
\frac{\partial A\overline{f}}{\partial t}=-\int_{\Gamma_{A(z)}}f(\mathbf{u_h}\cdot\mathbf{n})dl - \frac{\partial A\overline{w}\overline{f}}{\partial z} - \frac{\partial A\overline{w'f'}}{\partial z} - \frac{\partial A\overline{F_z}}{\partial z} + B_{F_z} + A\overline{R_f}.
\end{equation}
At this stage it is timely to distinguish between turbulent and non-turbulent flux, namely $F_z=F_{tz}+F_{nz}$, and define "effective" turbulent flux, $\overline{F_{tz}^*}=\overline{F_{tz}}+\overline{w'f'}$. This effective turbulent flux includes horizontally-averaged small-scale turbulent flux $(\overline{F_{tz}})$ and the flux mediated by organized flow structures, $\overline{w'f'}$. We also assume that the total non-advective flux $F_z$ at the bottom is the same at all bottom locations of the depth $z$. Then, taking into account the above hypotheses and

\begin{equation}
B_1=\int_{x_1(z)}^{x_2(z)}\left[ \frac{\partial y_2}{\partial z} - \frac{\partial y_1}{\partial z} \right]dx=\frac{dA}{dz},
\end{equation}
we transform \eqref{eq:fbasicav2} to

\begin{equation}
\label{eq:fbasicav3}
\frac{\partial A\overline{f}}{\partial t}=-\int_{\Gamma_{A(z)}}f(\mathbf{u_h}\cdot\mathbf{n})dl - \frac{\partial A\overline{w}\overline{f}}{\partial z} - \frac{\partial A\overline{F_{nz}}}{\partial z} - \frac{\partial A\overline{F_{tz}^*}}{\partial z} + \frac{dA}{dz}(F_{nz,b}(z)+F_{tz,b}(z)) + A\overline{R_f},
\end{equation}
where $F_{*,b}(z)$ denote bottom values of fluxes at depth $z$. The mean vertical velocity, $w$, may be expressed from the horizontally integrated continuity equation \eqref{eq:contin}:

\begin{equation}
\frac{\partial A\overline{w}}{\partial z}=B_w-\int_{\Gamma_{A(z)}}(\mathbf{u_h}\cdot\mathbf{n})dl,
\end{equation}
where $B_w \approx 0$ according to assumption of quasi-horizontal bottom. This means, $w$ arises from disbalance between inflows and outflows and subsequent water level change. For the LAKE model hasn't been applied for water bodies with significant water level change, the term with $w$ is omitted in \eqref{eq:fbasicav3} in the model equation set. The next strong assumptions applied to deliver a familiar 1D lake model equation from \eqref{eq:fbasicav3} are:

\begin{itemize}
  \item the non-turbulent flux, that is the shortwave radiation flux if $f$ is a temperature, is horizontally uniform, so that $\overline{F_{nz}(z)}=F_{nz,b}(z)=F_{nz}(z)$;
  \item the 'effective' turbulent flux may be represented via the gradient of mean quantity: $\overline{F_{tz}^*}=-k_{eff}\frac{\partial \overline{f}}{\partial z}$;
  \item the source averaged horizontally, $\overline{R_f(f,...)}$, may be approximated as the same function of mean values, $\overline{R_f(f,...)}=R_f(\overline{f},...)$.
\end{itemize}

Substituting these statements into \eqref{eq:fbasicav3}, we finally get:

\begin{equation}
\label{eq:fbasicav3}
\frac{\partial \overline{f}}{\partial t}=-\frac{1}{A}\int_{\Gamma_{A(z)}}f(\mathbf{u_h}\cdot\mathbf{n})dl + \frac{1}{A} \frac{\partial}{\partial z}\left( A k_{eff} \frac{\partial \overline{f}}{\partial z} \right) - \frac{1}{A}\frac{\partial A F_{nz}}{\partial z} + \frac{1}{A} \frac{dA}{dz}(F_{nz}(z)+F_{tz,b}(z)) + \overline{R_f}.
\end{equation}


\chapter{Liquid water equations for snow cover}
\label{app:snow}

In the model, liquid water content, $W$, is expressed as the ratio of liquid water mass to the total snow mass, $W=\rho_{w,sn}/\rho_{sn}$. To derive an equation for $W$, \eqref{eq:snowwat}, first note that

\begin{equation}
\label{eq:Weq}
W'=\frac{\rho_{w,sn}'}{\rho_{sn}}-\frac{W\rho_{sn}'}{\rho_{sn}},
\end{equation}
where the primes denote time derivative, for brevity. The local mass balance of liquid water is governed by gravitational infiltration, $\gamma_v,~[\gamma_v]=m/s$, and by the rate of freezing, $F_{fr}^*$:

\begin{equation}
\label{eq:rhowsn}
\frac{\partial \rho_{w,sn}}{\partial t} = - \frac{\gamma_v \rho_{w0}}{\partial z} - F_{fr}^*,
\end{equation}
where we introduced the reference liquid water density, $\rho_{w0}$, to convert the volumetric flux to mass flux. In turn, the snow density is determined by gravitational water flux, and by the snow compaction, $C_{sn}$:
	
\begin{equation}
\label{eq:rhosn}
\frac{\partial \rho_{sn}}{\partial t} = - \frac{\gamma_v \rho_{w0}}{\partial z} + C_{sn}.
\end{equation}
Inserting \eqref{eq:rhosn} and \eqref{eq:rhowsn} in \eqref{eq:Weq} yeilds \eqref{eq:snowwat}.


\end{document}
