Mechanistic Modeling of the Variability of Methane Emissions from an Artificial Reservoir

Abstract

The mechanistic model LAKE2.3 was tested for its capability to predict of methane (CH₄) emissions from reservoirs. Estimates of CH₄ emissions from the Mozhaysk reservoir (Moscow region) provided by the model showed good agreement with instrumental in situ observations for several parameters of the water ecosystem. The average CH₄ flux calculated by the model is 37.7 mgC-CH₄ m⁻² day⁻¹, while according to observations, it is 34.4 mgC-CH₄ m⁻² day⁻¹. Ebullition makes the largest contribution to the emissions from reservoirs (up to 95%) due to low methane solubility in water and the high oxidation rate of diffusive methane flux. During the heating period, an increase in methane emission is observed both in the model and empirical data, with a maximum before the onset of the autumn overturn. An effective parameter for calibrating the diffusive methane flux in the model is the potential rate of methane oxidation. For ebullition flux, it is the parameter q10 (an empirical parameter determining the relationship between methane generation and temperature) because methane production in bottom sediments is the most important. The results of this research can be used to develop mechanistic models and provide a necessary step toward regional and global simulations of lacustrine methane emission using LAKE2.3.

1. Introduction

The most important greenhouse gases in the atmosphere directly affected by anthropogenic emissions are carbon dioxide and methane. Methane is important since its global warming potential per molecule is 72 times higher than that of carbon dioxide over 20 years [1]. In addition, the relative growth rates of methane content in the atmosphere significantly exceed those for CO₂ and N₂O: an increase in the average global surface methane concentration relative to the pre-industrial period is 167% versus 46% for CO₂ and 24% for N₂O [2].

There are natural and anthropogenic sources of methane in the atmosphere. The most important natural sources include wetlands, tropical forest ecosystems, and lakes. Anthropogenic sources of methane in the atmosphere are landfills, agriculture, especially rice paddies and cattle pastures, industry, etc. [3]. Artificial reservoirs are also a notable anthropogenic source of methane for the atmosphere.

The main source of methane in the water column of natural lakes and reservoirs is the anaerobic decomposition of organic matter in sediments (Figure 1). Even a small presence of oxygen in the bottom horizons can inhibit the activity of archaea—the main agents of organic matter (OM) anaerobic destruction [4]. The dominating biogeochemical pathways of methane formation in sediments are hydrogenotrophic and acetoclastic reaction chains (Equations (1) and (2), respectively).

$$4{\mathrm{H}}_2^{+}+\mathrm{C}{\mathrm{O}}_2^{}=\mathrm{C}{\mathrm{H}}_4^{}+2{\mathrm{H}}_2\mathrm{O}$$

(1)

$$2\mathrm{C}{\mathrm{H}}_2\mathrm{C}\mathrm{O}{\mathrm{O}}^{-}+2{\mathrm{H}}^{+}=\mathrm{C}{\mathrm{H}}_4+\mathrm{C}{\mathrm{O}}_2+2\mathrm{H}\mathrm{C}{\mathrm{O}}_3^{-}$$

(2)

Typically, in the upper part of the bottom sediments, the acetoclastic chain of OM decomposition prevails, and in deeper layers, the hydrogenotrophic pathway provides the main contribution [5]. The composition of the OM settling at the surface of the bottom sediment is important—with the arrival of labile (rapidly decomposing) organic matter, methane production is significantly accelerated [6]. The distribution of methane fluxes and transformations in an artificial reservoir is described in detail in [7].

Methane produced in bottom sediments can reach the water–atmosphere interface in the aqueous form (diffusive flux of dissolved CH₄) or as a free gas in bubbles (ebullition process). In addition to the formation of methane directly in the bottom sediments, a significant contribution to the whole-lake CH₄ production can be provided by the decomposition of macrophytes in shallow waters [8]. The diffusive flux depends on the dissolved methane concentration gradient as well as the effective diffusion coefficient (including molecular and turbulent transport). When entering the oxygenated water layers, dissolved methane is oxidized by methanotrophic microorganisms, with a rate depending on concentrations of both methane and oxygen [9]. Up to 90% of the diffusive flux can be oxidized between the oxycline and the water surface [10]. Due to the composition of the above-mentioned factors, the CH₄ content increases from the surface to the bottom layers [11].

The key feature of ebullition flux is that it reaches the water surface much faster than the diffusive counterpart and cannot be directly oxidized. The bubble flux of methane depends on the local depth of the reservoir, as well as the water level dynamics. Artificial reservoirs are characterized by sharp level decreases during water flow regulation. A sharp drop in hydrostatic pressure at the bottom leads to an increase in bubble volume in sediments [12]. The local reservoir depth affects the magnitude of the bubble flux due to the dissolution of gas bubbles in the water column during their ascent, especially large-diameter bubbles [13]. With a smaller distance, the bubble passes from the bottom to the surface, and the fraction of gas molecules that have passed into the dissolved phase and subsequently oxidized in the oxycline is smaller [14].

An important factor affecting the total methane flux from the reservoir to the atmosphere is its trophic status. An increase in the phosphorus load to the reservoir, which leads to the eutrophication process, and related amplification of chlorophyll content in the water lead to an increase of methane emissions from reservoirs by 30–90% and may even bring lakes and reservoirs closer to swamps in terms of annual methane emissions into the atmosphere [15]. For eutrophic reservoirs, an essential seasonal increase in methane emissions takes place during the algae blooming period, when the rapid development of phytoplankton leads to a significant increase in the primary production of OM [16]. Currently, the detritus sedimentation rate and, hence, the amount of incoming labile OM significantly increase, leading to rapid depletion of oxygen in the bottom layer and supplying fresh substrate for anaerobic OM decomposition to methane.

In the inventory of methane emissions from reservoirs, it is important to take into account the horizontal components of the methane fluxes. A significant contribution to gross methane emissions from many reservoirs associated with hydroelectric power plants is provided by methane degassing in the water flow, leaving the water body through the dam into the downstream. For deep reservoirs with small residence time, the methane output during degassing can account for up to 70% of the total emission to the atmosphere [17].

According to available estimates, methane emissions from artificial reservoirs range from 2 to 122 Tg year⁻¹, or from 0.5 to 10% of the total methane flux into the atmosphere from the Earth’s surface [18,19,20,21,22,23,24]. Estimates of global methane emission from reservoirs vary significantly due to the differences in assessment methodologies, as well as in the data sets involved. The assessment methods typically assume extrapolation of statistical relationships obtained for the multitude of studied objects to the remaining (unexplored) reservoirs, using a set of predictors such as climatic zone, average depth, age, etc. [21]. This approach, being useful for global estimates, may fail in the prediction of methane emissions from particular reservoirs and even regions.

Current methane emission measurement technologies, such as “floating chambers”, eddy covariance, etc., have a few limitations: “floating chamber” measurements cannot be applied during unpleasant weather and need continuous series of field observations, eddy covariance has an uncertainty of results due to landscape footprint problem. Comparison and more detailed descriptions of different observation methods are shown in [18]. A more reasonable methodology for estimating methane emissions from the surface of water bodies not covered by observations or covered by a limited set of measurements involves mechanistic modeling of the key processes of production, consumption, transport, and emission of methane from a reservoir. The model needs to be tested and calibrated using the available empirical data on the spatial and temporal heterogeneity of methane fluxes in the water body of interest. After determining the optimal values of the model parameters and simulating the time series of methane fluxes at the water—atmosphere interface, it is possible to obtain refined estimates (inventory) of the annual methane flux for reservoirs, given their main morphometric, hydrological characteristics and atmospheric forcing.

One of the process-based models capable of performing such estimates is the one-dimensional thermohydrodynamic and biogeochemical model LAKE [25,26]. Models of this type have been previously used to estimate CH₄ and CO₂ emissions from natural lakes with high residence time [27,28,29,30,31,32,33]. However, artificial water bodies are often characterized by significant horizontal heterogeneity in the distribution of both physical and biogeochemical variables, presenting a challenge for 1D models to successfully reproduce surface energy and mass exchange [34]. Therefore, one of the objectives of this study is to assess the applicability of a 1D (vertical) approach to simulate the concentration and fluxes of methane in reservoirs using the LAKE model.

In this study, the Mozhaysk artificial reservoir (Moscow region, Russia) was chosen as an object to test the capability of the LAKE model to simulate methane fluxes and to apply the concept of improved inventory described above. The paper compares in situ observed CH₄ concentrations and fluxes (as well as water temperature and dissolved oxygen) to results of the LAKE model and discusses possible ways to improve the quality of modeling and thus arrive at more accurate estimates of methane emissions from artificial reservoirs.

This paper includes an Introduction, providing a short overview of the processes responsible for methane fluxes in the water column and on the surface of the freshwater body; Materials and methods, describing the Mozhaysk reservoir, observation methods and results, as well as the LAKE model version 3, adopted for reservoirs with high throughflow and surface level variations; Results and discussion, presenting results of water temperature, oxygen and methane emission simulation by LAKE 3.0, comparison of simulations with empirical data, and quantification of the influence of selected factors on methane emission based on sensitivity experiments with the model.

2. Materials and Methods

2.1. Study Site

The Mozhaysk artificial reservoir was used in this study as the site of field research and verification object for the LAKE model. This is a small morphologically simple (a reservoir without flooded valleys of large tributaries) valley reservoir in the Moscow region (55.5948° N, 35.8221° E) with a slow water exchange (

Table 1. Morphological characteristics of the Mozhaysk reservoir (all characteristics are given for the mean headwater level) [35].

Length, km	Max Width, km	Mean Width, km	Max Depth, m	Area, km2	Volume, km3	Range of Level Fluctuations, m/Year	Water Residence Time, Year
28	2.6	1.1	22.6	30.7	0.235	6	0.6

).

To study the spatial heterogeneity of methane fluxes in the reservoir, five reference observation stations were selected over the flooded riverbed (Figure 1). This arrangement of measurement stations at a uniform distance from each other in the longitudinal direction of the reservoir facilitates the study of physical and biogeochemical characteristics distribution from upstream (affected by the tributaries) through the transformation zone in the middle course of the reservoir and to the lower part with a calm lake regime. In addition to observations over the flooded riverbed of the Moscow River, measurements were also carried out on other morphological sites within the selected reservoir compartments I–V (Figure 1), specifically over the flooded floodplain of the river valley.

2.2. Field Observations

The methane flux was measured by the floating chamber method [7,36]. The method is based on measuring the concentration difference in the hermetic chamber, which is installed on the water at the beginning and the end of the exposure period. We used chambers of two types: a common chamber determining the sum of the bubble and diffusion flux and a diffusion chamber, which is equipped with a screen that cuts off the bubble flux. The time of chamber exposition varied from 30 min to 1.5 h. Water samples for detecting dissolved methane using the “headspace” method [7] were taken simultaneously with the measurement of methane fluxes by chambers. Samples were taken from the surface and bottom water horizons, as well as above and below the layer of the largest density gradient in the presence of stable temperature stratification. Sampling was accompanied by measurements of the water temperature, conductivity, and dissolved oxygen with YSI ProODO and Pro30 probes (YSI Inc., a Xylem brand, Yellow Springs, OH, USA). In addition, the air temperature, atmospheric pressure, and wind speed were measured by Davis Vantage Pro meteorological suit at Station IV from May to September. Automatic measurements of water temperature were carried out at Station IV by temperature loggers at depths 0.5, 1, 2, 3, 4, 5, 7, 10, and 14 m every 15 min (Hobo pendant MX 2202 and Hobo Water Temperature Pro U22-001 (Onset Comp., Bourne, MA, USA)). Dissolved oxygen can be measured by HOBO logger every 3 h on the surface layer (0.5 depth) and bottom layer (14 m). For vertical measurements, the YSI ProODO with a 30 m cable was also used.

Measurements were divided into two types—measurements at all 5 stations (Figure 1) (carried out for one day) and more frequently temporal measurements on Station IV (carried out on different time intervals, in some dates 2–4 times during the day) for more detailed information about temporal variability of methane fluxes and concentrations.

Table 2. Dates of methane concentrations and fluxes measurements 2015–2020.

Year	Measurements at All 5 Stations	Additional Flux Measurements at Station IV
2015	-	25 June; 29 June; 8–10 July; 28 July
2016	24 June; 4 July; 13 Jule; 22 August Only concentrations	26 June; 27 June; 3 July; 9 July; 12 July; 27 July; 28 July; 21 August
2017	4 July; 5 July; 20 August Only concentrations	20 June; 24 June; 27 June; 28 June; 3 July; 10 July; 28 July; 29 July; 31 July; 5 August; 19 August; 2 September; 9 September
2018	24 June; 7 July; 19 August	20 April; 9 May; 25 May; 17 June; 19 June; 26–29 June; 2 July; 8–10 July; 28 July; 17 August; 22 August; 28 August; 2 September; 8 September; 16 September; 21 September
2019	10 June; 24 June; 5 July	22 June; 23 June; 26 June; 2 July; 7 July; 9 July; 12 July; 12 August; 9 August
2020	25 May; 12 August	16 June; 24 June; 9 July; 12 July

shows the dates of measurements. Measurements were predominantly carried out during the summer period of straight stratification.

The errors of flux and concentration values were evaluated to obtain more accurate measurement results. Relative error was calculated as a result of equivalent measurement series. During these measurements, 3–12 successive samples of water methane concentration from the same location and depth and flux from the same chamber with 20-min intervals were taken. The results of this experiment showed that the relative error has no relationship with the value of the flux or concentration and no relationship with replication numbers. Therefore, the maximum value of relative errors was applied as a method of uncertainty estimation (

Table 3. Experiment of relative error calculation for methane concentration and methane flux measurement (bold shows the measurements with the highest relative error).

Type of Measurements	Number of Replications	Mean Value	Relative Error, %
Concentration	3	11.3 μLCH4 L−1	13.9
Concentration	6	578.3 μLCH4 L−1	15.3
Concentration	8	18.7 μLCH4 L−1	1.8
Concentration	4	27.7 μLCH4 L−1	8.2
Concentration	5	18.0 μLCH4 L−1	2.2
Concentration	7	1165.5 μLCH4 L−1	13.9
Flux into atmosphere	7	15.2 mgC-CH4 m−2 d−1	27.0
Flux into atmosphere	7	3.09 mgC-CH4 m−2 d−1	25.0
Flux into atmosphere	12	6.74 mgC-CH4 m−2 d−1	25.7

).

The relative error values for measurements are 15.3% for methane concentration and 27.0% for methane flux.

2.3. Satellite Measurements of Water Surface Temperature

With the advent of medium- and high-resolution satellites, it has become possible to detect various characteristics of any size inland water bodies [37]. A multispectral method using data from several satellite missions was used to determine the remote surface temperature of the Mozhaysk reservoir. When constructing the reservoir mask, only those pixels that are reliably located on the water surface and do not capture areas of the shoreline were considered [38], and cloudy pixels were also excluded. The methodology for cloud pixel extraction is described in [39]. Water temperature data from Landsat-8 satellite data were obtained using a two-channel algorithm [40]. This algorithm allows the reconstruction of the real surface water temperature. The Landsat-8 satellite carries the TIRS (Thermal InfraRed Sensor) instrument, which measures in two channels in the far infrared range (10.30–11.30 µm; 11.50–12.50 µm), which allows the provision of atmospheric correction of measurements. The signal received by the satellite is the sum of signals from the water surface (temperature proper) and from the atmosphere (related to atmospheric transmittance). The emission characteristics of the atmosphere in the two channels are different, while the emission characteristics of the surface are the same (in the first approximation). Thus, the water surface temperature ($T_w$), according to Landsat-8 data, can be calculated according to the following formula:

$$T_w=1.911T_{10}-0.8554T_{11}+0.7837$$

(3)

where $T_{10}$ and $T_{11}$ are temperatures obtained from channels 10 and 11, respectively. Calibrated temperature data from MODIS (Moderate Resolution Imaging Spectroradiometer) measurements were used to determine the coefficients. Sentinel-1 (2016) and Sentinel-3/SLSTR (2017) data were also used. The lake water temperature was calculated using the algorithm presented in [41,42]. The surface temperature using the MODIS scanning spectroradiometer on the TERRA and AQUA satellites was determined using the methodology outlined in [42]. The L1C product of the MIRAS microwave radiometer (SMOS/ESA satellite) was also used to obtain the seasonal dynamics of the brightness temperature of the cell corresponding to the reservoir. Periods of characteristic open water and ice cover brightness temperature values were identified [43].

2.4. The Summary of the LAKE Model

The LAKE is a one-dimensional (with a partial representation of horizontal inhomogeneity) thermohydrodynamic model of a water body with a module of biogeochemical processes responsible for the formation of methane and carbon dioxide (https://mathmod.org/lake/, accessed on 20 October 2023). Relevant model code and all scripts can be found on the “zenodo” database [44]. The model has been used previously to simulate the thermodynamic regime of reservoirs, as well as the CO₂ [45] and CH₄ [26] dynamics. A version of the model LAKE2.3, summarized below, was used to estimate the CH₄ flux from the reservoir.

The basic assumption of the model is that the system of equations for horizontally averaged physical and biogeochemical quantities of the water can be closed by reasonable mechanistic hypotheses, containing parameters enabling efficient model optimization in terms of target variables. In physical terms, it means that the vertical profiles of prognostic variables and total fluxes to the atmosphere can be successfully represented without explicitly simulating the lateral heterogeneities imposed by 3D circulation. The target variables include basic thermodynamic characteristics and greenhouse gases (concentrations and fluxes).

This section provides a general description of the mathematical model with links to publications where individual modules are given in more detail. Particular attention is paid to the representation of physical mechanisms and biogeochemical effects associated with the water flow through the reservoir since their representation mainly distinguishes the 2.3 version of the model from the preceding ones.

The one-dimensional (in vertical) equations for thermodynamic, dynamic, and biogeochemical variables in a reservoir are the result of an averaging operation over the horizontal section of a water body (Figure 2) applied to three-dimensional equations for horizontal momentum components, the continuity equation, the heat equation, the dissolved/solid species balance equations.

The resulting form of a one-dimensional equation for a scalar value f (including any of the horizontal velocity components u, v) in an incompressible fluid has the form:

$$\frac{\partial \overline{f}}{\partial t}=-\frac{1}{A}{\oint }_{{\Gamma }_A}^{}f\mathit{u}·\mathit{n} dl-\frac{1}{A}\frac{\partial A\overline{w}\overline{f}}{\partial z}+\frac{1}{A}\frac{\partial }{\partial z}\left[A\left(k_T+k_m\right)\frac{\partial \overline{f}}{\partial z}\right]-\frac{1}{A}\frac{\partial A{\overline{\Phi }}_f}{\partial z}+\frac{1}{A}\frac{dA}{dz}\left(F_{f,{\Gamma }_{A\left(z\right)}}+{\Phi }_{f,{\Gamma }_{A\left(z\right)}}\right)+\overline{R}$$

(4)

where z—the vertical coordinate directed along gravity, with the origin at the water surface, t—time, A(z)—the horizontal section, ${\Gamma }_{A\left(z\right)}$—the closed boundary of the section A(z), dl—the element of the ${\Gamma }_{A\left(z\right)}$ boundary length (Figure 2), n—the external normal to ${\Gamma }_{A\left(z\right)}$, u = (u, v)—the horizontal component of velocity, w—the vertical projection of velocity, $F_f$—the total diffusion flux of f due to turbulence and molecular exchange, ${\Phi }_f$—the sum of non-diffusion and non-advective fluxes of f (for example, radiation flux in the equation for temperature, bubble flux in the equations for dissolved gases),$F_{f,{\Gamma }_{A\left(z\right)}},{\Phi }_{f,{\Gamma }_{A\left(z\right)}}$—values of the corresponding fluxes at ${\Gamma }_{A\left(z\right)}$, i.e., at the bottom of depth z (within ${\Gamma }_{A\left(z\right)}$ these fluxes are assumed to be constant), $k_T$ and $k_m$ are the coefficients of turbulence and molecular diffusion, R—the sum of all terms of the original three-dimensional equation, except for the total derivative and divergence of fluxes (i.e., sources and sinks in the equations for biogeochemical substances, pressure gradient, and Coriolis force in the equations of motion, etc.), overbar $\overline{\left(\dots \right)}$ stands for horizontally averaged expression.

Equation (4) is valid for the case of a bottom with small slopes (so that horizontal components of diffusion fluxes and vertical velocity at the bottom are negligible), as well as for a reservoir formed by vertical walls with zero normal diffusion fluxes and a horizontal bottom. In the context of this study of a water flow through the reservoir, the first and second terms on the right-hand side of (4) are of particular interest since they are responsible for the inflow by the tributaries and removal by the discharge of the value f and advection of f by the average vertical velocity, respectively.

The equations in the form (4) are solved in the model for the following variables:

−

horizontal speed components;

−

temperature;

−

salinity;

−

concentration of dissolved oxygen and methane;

−

concentration of carbon in the following forms: living organic particles (phyto- and zooplankton), dead organic particles (detritus), autochthonous and allochthonous dissolved organic carbon, dissolved inorganic carbon;

−

concentration of phosphorus of the dissolved inorganic phosphorus (phosphates).

Equations of the form (4) are supplemented by boundary conditions at the upper and lower boundaries. For temperature, the heat balance equation is used at the upper boundary (z = 0), and the continuity of temperature and heat flux is at z = H (H is the maximum depth of the reservoir). The heat balance on the surface is calculated under a given time series of shortwave and longwave radiation fluxes and basic meteorological parameters; the turbulent fluxes of sensible, latent heat and momentum are calculated using the Monin-Obukhov similarity theory [46,47].

The equation for the mean vertical velocity is obtained by averaging over the horizontal cross-section of the continuity equation:

$$\frac{\partial A\overline{w}}{\partial z}=-{\oint }_{{\Gamma }_A}^{}\mathit{u}·\mathit{n} dl,$$

(5)

given the impermeability condition ${\left.\overline{w}\right|}_{z=H}=0$.

Integrating (4) and taking into account water fluxes at the surface, the equation for H is obtained:

$$\frac{dH}{dt}=r-E-\frac{1}{A\left(0\right)}{\int }_0^H{\oint }_{{\Gamma }_A}^{}\mathit{u}·\mathit{n} dldz+M,$$

(6)

which expresses the change in the water level as a result of the imbalance of the discharges of inlets and outlets, precipitation r, evaporation E from the surface, as well as the processes of freezing and melting of ice and snow cover (represented by the summand M).

The system, which consists of one-dimensional equations of the form (4) and the continuity Equation (5), is closed using additional hypotheses and parameterizations [25]. Thus, the turbulent closure k-ε is used to calculate the coefficients of turbulent viscosity and thermal conductivity (diffusivity); the expression for the coefficient of thermal conductivity (diffusivity) also includes the additive coefficient of “background diffusion”, representing vertical mixing due to the internal waves breaking and other mixing effects not accounted for in standard turbulent closures. Parameterizations of the fluxes of scalar quantities and momentum at the bottom surface $F_{f,{\Gamma }_{A\left(z\right)}},{\Phi }_{f,{\Gamma }_{A\left(z\right)}}$ at z < H also play an important role. Heat and methane fluxes at the bottom are calculated by solving one-dimensional problems for temperature and methane in layers (columns) of bottom sediments with a boundary with water bodies at different depths [25], including the maximum depth, so that the sum of these boundaries is the entire bottom area of the reservoir (Figure 2). Momentum flux through the bottom surface $\left(F_{u,{\Gamma }_{A\left(z\right)}},F_{v,{\Gamma }_{A\left(z\right)}}\right)$ is calculated according to the linear or quadratic law with respect to average velocity $\left(\overline{u},\overline{v}\right)$ with a calibration multiplier.

The one-dimensional model of heat and moisture transfer in bottom sediments takes into account the possibility of phase transitions of water. The equation for methane concentration in sediments includes production, molecular diffusion, and removal of methane in the form of bubbles when the critical CH₄ concentration is exceeded [27].

Non-diffusive and non-advective fluxes ${\Phi }_f$—are the kinematic shortwave radiation flux ${\Phi }_T=S$ in the heat equation (computed according to the Beer-Buger-Lambert law for the infrared, photosynthetically active, near and far infrared wavelength ranges) and the bubble flux in the equations for the concentrations of dissolved oxygen, carbon dioxide, and methane. The bubble model is based on parameterizations from [48] and takes into account the exchange of five gases between the water and the bubble: O₂, CO₂, CH₄, N₂, and Ar.

A separate task is to construct parametrizations for the terms $\overline{R}$ in one-dimensional models. Thus, parameterization of the averaged horizontal pressure gradient in the momentum equations in the LAKE model allows the reproduction of seiches with horizontal wave number 1 [26]. Sources and sinks $\overline{R}$ in the equations for biogeochemical substances are set using parameterizations from [29,49,50,51,52,53] and take into account the following processes:

−

photosynthesis, respiration, exudation, and death of phyto- and zooplankton;

−

aerobic decomposition of dissolved organic compounds and detritus;

−

photochemical decomposition of dissolved organic compounds;

−

aerobic oxidation of methane.

The consumption of dissolved oxygen by bottom sediments with the release of dissolved inorganic carbon and dissolved inorganic phosphorus are expressed by terms $F_{f,{\Gamma }_{A\left(z\right)}}$ in the balance equations of these variables according to [54].

The one-dimensional model with a partial representation of horizontal inhomogeneity presented above allows the calculation of the following fluxes of gases into the atmosphere:

−

diffusion flux from the surface of the reservoir, where the gas exchange coefficient is provided according to the surface renewal model [55];

$$F_C=k_{ge}\left({\left.C\right|}_{z=0}-C_{ae}\right),$$

(7)

where $k_{ge}$ is the gas exchange coefficient (piston velocity), m/s, ${\left.C\right|}_{z=0}-C_{ae}$ is the difference between air and surface water gas concentrations;

−

bubble flux of methane, carbon dioxide, and oxygen are calculated separately above bottom sediments located at different depths (Figure 2);

$$F_{B,f,k}=M_{f,k}\times n_{b,k}\times v_{b,k},f={CH}_4,{CO}_2,O_2,$$

(8)

where $M_{f,k}$ is a molar gas concentration, $n_{b,k}$—bubbles concentration in water, $v_{b,k}$—vertical velocity of bubbles, $k$ is the index of sediments column;

−

advective flux through the outflow (or turbines in the case of hydroelectric power plants) (Equation (9)).

$$F_{degas}=Q\times (C_{upstream}-C_{downstream}),$$

(9)

where $Q$—is water discharge through turbines, $C_{upstream/downstream}$—methane concentration in a reservoir and below a dam, respectively.

The model includes multilayer modules for calculating the transfer of heat and liquid water in layers of ice and snow [56].

The system of the model equations is solved by finite-difference methods. The model uses a Crank–Nicolson center-differences scheme for diffusion terms, an explicit Euler scheme for sources/sinks of concentrations, and a MUSCL scheme for vertical advection [57].

2.5. Setup of Numerical Experiments

The input data to the model is divided into hydrological information (discharge of water, heat, momentum, and water constituents by main tributaries, water level or discharge through a dam) and meteorological forcing (temperature, humidity, atmospheric pressure, wind speed, shortwave, and longwave radiation, precipitation). Water level data were extracted from the Mozhaysk Hydroelectric Power station archive, and ERA5-Land reanalysis was used as a source of meteorological data. These meteorological data were corrected according to the available series of ground observations both at the reservoir and at the nearest meteorological stations, Mozhaysk and Naro–Fominsk [58]. Data were corrected for the modeling period of 2015–2019.

The baseline numerical experiment with the LAKE model was carried out for the period from 2015 to 2019, whereas the period November–December 2014 was used as the model spin-up period. The time step of the model finite-difference scheme was set to 10 s, the step of meteorological data was 1 h, and the input hydrological information was updated daily. The water temperature observed in the fall of 2014 was taken as the initial conditions. The oxygen saturation degree at the initial moment was set to 100%, which is realistic for the autumn well-mixed homothermy stage. The initial methane concentration was set to 0 since no corresponding observations were available at that time. A total of 22 horizontal levels were set in the model with an approximate spacing of 1 m, which corresponds to the maximum depth of the reservoir at the initial instant, 5 columns of bottom sediments evenly distributed over the depth of the reservoir, each 1 m thick, with the grid step inside each column 10 cm.

Model parameters were measured in situ where possible or set up according to literature data. For example, the extinction coefficient of solar radiation was set as constant according to Secchi disk measurements on the reservoir, which equals 1.4 m⁻¹. The coefficient of additional diffusion in the thermocline is important for the correct simulation of the diffusive heat flux in metalimnion with stratification strong enough to cause turbulent moments vanishing in standard turbulence closured (including k-ε closure used in the LAKE model). This coefficient is semi-empirical and is expressed by the formula:

$$k_{BD}=k_0{A}_0^{k_1}{\left(N^2\right)}^{{-k}_2},$$

(10)

where $k_0$ = 8.17 × 10⁻⁴, $k_1$ = 0.56, $k_2$ = −0.43 are empirical coefficients, $A_0$ is the water surface area, and N is the Brunt–Väisälä frequency [25,59].

The model output step is 1 h. All the variables computed on the model temporal grid spacing 10 s were averaged to hourly intervals. The hourly values were then averaged over monthly and annual intervals. The LAKE model inputs—temperature, wind speed, pressure, and humidity—were linearly interpolated to 10 s mesh, precipitation, and radiation fluxes were interpolated piecewise constant; the river inflow data were treated piecewise constant as well.

3. Results and Discussion

3.1. Methane Emission by In Situ Measurements

The greatest temporal variability of methane fluxes from the Mozhaysk reservoir is typical for the summer season. Results of measured methane fluxes during the summer period are presented (Figure 3). Information on this variability was collected during observations conducted at Station IV (Figure 1) for the period 2015–2019.

For most years, there is a similar tendency for the methane flux to increase during summer (2016–2018). The maximum flux value is typically reached at the end of the summer period, i.e., at the commencement of the autumn convection. By this time, a thick anaerobic zone with dissolved O₂ content of less than 1 mg L⁻¹ is formed in the hypolimnion due to the stable temperature stratification [60]. This favors methane accumulation in the bottom layers of water. In addition, by the end of summer, the productivity of the reservoir significantly increases, and a large amount of dead organic matter (detritus) subsides to the bottom sediments, thus adding a substrate for anaerobic destruction (in agreement with the established link of trophic status to methane emissions [15]). This, in turn, contributes to a significant increase in the total methane flux. At high values of the total methane flux, the ebullition flux can account for more than 90% of the total emission [60]. It is worth noting that the increase of CH₄ flux during summer may occur gradually (as in 2017, when the temperature stratification was weak) or include a sharp release of methane when a stable temperature stratification took place (as at the end of summer 2018).

The first step of methane emission estimation was to estimate the flux from the riverbed morphometric part of the reservoir. Quasi-synchronous (conducted for 1 day) observations on the flooded riverbed reservoir part of total methane flux in different parts of the Mozhaysk reservoir (stations I–V) (Figure 1) were carried out in the summer periods of 2018–2020. The average ratio of CH₄ flux at riverbed Station IV to the zonally averaged flux from the entire riverbed reservoir area was calculated based on the results of all the surveys (Figure 4a).

For horizontal averaging of CH₄ riverbed flux, the areas of the reservoir compartments (Figure 1) were used. For each compartment, a methane riverbed flux value was assigned from the corresponding observation station (Figure 1). Then, the following regression was assumed between methane riverbed flux averaged on all Mozhaysk reservoirs and methane flux on Station IV. Thus, the multiplier 1.3 was used to estimate the total emission from the entire riverbed area of the reservoir by the values of CH₄ flux obtained at Station IV only, where the flux on flooded river channel was measured more often (61 measurements at Station IV versus 8 whole-reservoir surveys). However, in Figure 4a, one value (August 2019) is significantly out of the general pattern. This is due to the very low water level of the reservoir in August 2019, at which the spatial distribution of methane fluxes differs greatly from the distribution observed during the typical water levels; therefore, at water levels below 179 m (Baltic Elevation System), the mentioned conversion multiplier from value at Station IV to entire riverbed reservoir area was assumed to be 0.4.

The values of methane fluxes considered above refer to the flooded river channel of the Mozhaysk reservoir since all main measurement stations (Figure 1) were in the deepest points of sections. The second step of emission calculation is to relate the values obtained for the riverbed part with the floodplain part of the Mozhaysk reservoir. Several quasi-synchronous measurements were carried out in each compartment and at different periods above the riverbed and above the floodplain in the same reservoir section (Figure 4b) to consider inundated floodplain areas, which have smaller depths and larger areas. The ratio of methane fluxes over the riverbed and floodplain differs significantly in all sections and between periods. There is no significant relationship between riverbed and floodplain methane fluxes. However, because these synchronous measurements covered the entire water area of the Mozhaysk reservoir, the average ratio of the floodplain flux to the riverbed flux for all measurements was taken as a conversion factor, which was assumed to be 0.38.

Thus, based on the temporal variability of methane fluxes at Station IV using the quantitative relationships of fluxes in different parts of the reservoir described above, we computed the time series of zonally averaged values of the total methane flux from the surface of the Mozhaysk reservoir based on measurements at Station IV. These time series were linearly interpolated between measurements at Station IV to obtain a series with daily resolution. The averaged fluxes of the daily resolution were multiplied by the area of the Mozhaysk reservoir for each day, thus obtaining the value of the daily methane emission. The sum of daily emission values provided the annual methane emission from the surface of the Mozhaysk reservoir.

The emission of methane during degassing downstream of the hydroelectric powerplant was estimated to be degassing through turbines. It was calculated as the difference in methane concentrations in the horizon of water intake in the upstream and downstream multiplied by the turbine discharge [17]. Thus, the values of total methane emissions during summertime for 2017–2019 (less measurements data of the methane flux during the summer period of 2015 and 2016 do not allow the estimation of emission in these years) from the Mozhaysk reservoir were obtained (

Table 4. Estimates of the annual methane emission from the Mozhaysk reservoir for 2017–2019 according to field observations.

Year	Emission from the Water Surface, tC-CH4	Degassing through Turbines and Downstream, tC-CH4	Total Emission, tC-CH4
2017	334	13	347
2018	256	9	265
2019	378	11	389

).

During the ice period, methane emission is assumed zero; in spring and autumn, when the reservoir is characterized by homogeneity of temperature and dissolved/solid substances, low water temperature values, and low productivity, the methane flux is set to a value 2.4 mgC-CH₄ m⁻² day⁻¹, which is the mean value of several observations during these periods in different years (September 2017, April 2018).

3.2. Thermodynamic, Oxygen, Methane Concentration Regime and Methane Emissions According to the LAKE Model

The quality of the model simulations was assessed in terms of 4 variables—water temperature, dissolved oxygen, dissolved methane, and methane fluxes at the water–atmosphere interface. The results of the LAKE model were compared with observations conducted during the summers of the 2016–2019 period at Station IV, neglecting the total emission estimations, which were compared to already averaged values (see Section 3.1).

3.2.1. Temperature Regime

The simulated water temperature time series were compared with the logger data from depths 0.5, 1, 2, 3, 4, 5, 7, 10, and 14 m at Station IV. The model reproduces the temperature variability in the upper mixed layer well enough, both at diurnal and seasonal scales (Figure 5).

Table 5. Statistical criteria for comparison of water temperature between model calculations and in situ measurements.

Surface Layer				10 m Layer
Year	Mean Delta, °C	Pearson R	RMSE, °C	Year	Mean Delta, °C	Pearson R	RMSE, °C
2016	−0.54	0.97	0.99	2016	1.32	0.91	1.89
2017	−0.84	0.96	1.20	2017	−0.51	0.96	0.98
2018	−0.01	0.92	1.05	2018	4.26	0.64	5.07
2019	0.00	0.98	0.97	2019	2.07	0.78	2.95
Average	−0.34	0.96	1.06	Average	1.77	0.66	3.09

shows the difference in mean values between calculated by model and measured data and the statistical criteria of data convergence for model temperature calculation quality. Mean delta, Pearson R, and RMSE were used as statistical criteria.

The temperature at 10 m depth is reproduced by the model less successfully and strongly depends on the conditions of a particular year. The best results were obtained for 2017 and 2016, while in 2018 and 2019, the model provided more intense seasonal vertical heat transport to deep layers compared to observed temperatures. The vertical heat transport in the metalimnion in the model is substantially governed by the coefficient of background thermal conductivity of a simplified form (Equation (4)), which does not allow mixing in metalimnion to be reproduced equally well under atmospheric forcing of different years. Figure 5 also demonstrates high-frequency oscillations of temperature at 10 m in the observed series, presumably originating from internal waves, not reproduced in the 1D model, as seiche oscillation parameterization [25] was switched off in this study.

The LAKE model reproduced the spatially averaged water layer temperature due to its 1D scale. Therefore, the model surface layer calculations were also compared with satellite observations described in Section 2.3 of this article. Satellite data allow the obtaining of averaged temperature data that can be averaged over the reservoir surface in contrast to station measurements. Water surface temperature obtained by the LAKE model and satellite measurements for 2016–2017 are visualized in Figure 6.

Model water surface temperature has even better data convergence with satellite observations than with in situ measurements. In general, model results of water surface temperature calculation are lower than the satellite data for 0.67 degrees for the 2016–2017 period. Statistical criteria are shown in

Table 6. Difference in mean values and statistical criteria for comparison of water surface temperature between model calculations and satellite observations.

Year	Mean Delta, °C	Pearson R	RMSE, °C
2016	−0.48	0.97	1.61
2017	−0.84	0.98	2.05
Average	−0.67	0.97	1.84

.

3.2.2. Dissolved Oxygen Regime

During summer, the measured dissolved oxygen concentrations in the surface water layer were significantly higher than those provided by the model (Figure 7a). In the LAKE model, the daily-averaged parameterization scheme of the carbon cycle processes is adapted from [50]. Therefore, the daily cycle of O₂ content in the surface layer is reproduced with suppressed amplitude. In the model, oxygen concentrations are in close equilibrium with the atmospheric O₂ content, which causes underestimation because water can be oversaturated with oxygen as a result of the photosynthesis process, according to field measurements. The smallness of deviations from the equilibrium state is explained by the approximate equality of sources (photosynthesis) and sinks (respiration, decomposition of dead organic matter) in the modeled epilimnion. In turn, this may be caused by the high rate of detritus oxidation, which does not leave time for detritus to sediment from the mixed layer due to gravitational deposition. However, significant model error in the O₂ concentration does not matter for the simulated methane oxidation rate in the surface layer. Oxygen concentrations of ~6–14 mg L⁻¹ are too high for limitation of the methanotrophic bacteria activity. In other words, both in the model and the observations, the oxygen content in the mixed layer significantly exceeds the half-saturation constant 0.33 mg L⁻¹ for O₂ in the Michaelis–Menten kinetics of CH₄ oxidation, used in LAKE [61].

For the processes of methane accumulation in the hypolimnion, which is important for CH₄ flux formation, a crucial variable is the oxygen content in the bottom layer, especially during the formation of an anoxic zone. The model satisfactorily reproduces the observed O₂ content at 14 m depth at Station IV for most periods of measurements (excluding the August 2017 oxygen drawdown), where such measurement data are available (Figure 7b).

The statistical criteria (Pearson R and RMSE) were not calculated for comparison between regular in situ measurements on Station IV and model results due to the irregular temporal frequency of these measurements. The only difference in mean was calculated for the period 2016–2019. For the surface layer, the average value of in situ measurements is 10.08 mg L^−1, and the mean of calculated values is 8.40 mg L⁻¹. Model concentrations were lower by 2.08 mg L⁻¹. For the bottom water layer, the average value of in situ measurements is 0.39 mg L⁻¹, and the mean of calculated values is 1.26 mg L⁻¹ (the difference in mean is 0.87 mg L⁻¹).

In parallel with single measurements, the oxygen concentration in water was obtained by logger gauges, described in Section 2.2. The comparison between model results and oxygen loggers is shown in Figure 8.

The model underestimated observed surface oxygen concentrations, but there is a strong agreement between modeled bottom water temperature and observed data (

Table 7. Difference in mean values and statistical criteria for comparison of oxygen content between model calculations and observations.

Year	Mean Delta, mgO2 L−1	Pearson R	RMSE, mgO2 L−1
2016 (surface)	−1.67	−0.26	2.71
2017 (surface)	−0.30	−0.03	2.09
2019 (bottom)	−0.34	0.99	0.52

).

3.2.3. Methane Concentration in Water

The LAKE model satisfactorily reproduces the observed temporal variability of the mixed-layer CH₄ concentration but overestimates the measured values on average by 0.3–0.5 μmol L⁻¹ (Figure 9). This can be partially related to the spatial variability of methane concentration. Observations (measurements at 5 stations described in Section 2.2) showed that methane surface concentrations could vary by 0.2–0.3 μmol L⁻¹ between the upper and lower parts of the reservoir (Stations I and V in Figure 1). Compared to CH₄ observations at the bottom layers at Station IV, the LAKE model demonstrates the lower rate of methane accumulation in hypolimnion (Figure 9b). However, seasonal concentration maxima coincide at the same time according to observations and the LAKE model. For comparison, Figure 9b shows the time series of methane concentration not only at a depth of 14 m (local depth of Station IV) but also at the lower model layer of 22 m (maximal depth of the reservoir). Even at the deepest horizon in the model, methane is accumulating less than what was observed at the local depth of Station IV. The reason for this discrepancy between the model and observations is likely to be two-fold. First, the 1D lake model by construction simulates an average CH₄ concentration, which should be lower than the methane content at the bottom, where the gas has its source. Second, the suggested link between mixed-layer productivity and methane production in sediments through sedimentation of detritus (see Section 3.1) is not represented in the current version of the LAKE model, which can also lead to some bias.

3.2.4. Methane Fluxes into the Atmosphere

The time series of the diffusion and ebullition methane fluxes averaged over the surface of the Mozhaysk reservoir simulated by LAKE are given in Figure 10. The main contributor of methane emission is ebullition—during the periods of maximum fluxes, it can contribute up to 95% of total emission into the atmosphere. The second component of the flux is diffusive emission, and the lowest values were attributed to methane degassing from turbines. During the winter period, the values of all 3 flux types are low. Values of more than 10 mgCH₄ m⁻² d⁻¹ can occur when the ice cover weakens. High emissions can occur in early spring, when the ice cover melts totally (up to 302 mgCH₄ m⁻² d⁻¹ 1 April 2018 and up to 405 mgCH₄ m⁻² d⁻¹ 28 February 2019). During the spring–summer warm-up period, methane fluxes increased slowly with high peaks of emission, which corresponds to reservoir level changes. The highest annual ebullition and diffusive fluxes were registered at the beginning of the autumn convective mixing stage when the density stratification in the water column weakens and disappears. Reservoir stratification break events can occur during the warm summer period because of strong wind–wave impact or also with water level drawdown events (high peaks of the flux appear during water level minimum extremes). Methane fluxes and water levels are shown in Figure 10.

In 2018 and 2019, the beginning of the autumn mixing stage and water level drawdowns happened at the same time, and the highest values of the flux of all calculation periods were registered (649 mgCH₄ m⁻² d⁻¹ 17 October 2018 and 1903 mgCH₄ m⁻² d⁻¹ 17 September 2019). The high temporal variability of the diffusive flux corresponds to methane surface concentration variations. Ebullition variations correspond to atmosphere pressure changes.

For comparison to the measured total CH₄ flux at Station IV, we used the sum of the bubble flux emitted from the second-deepest sediments column in the model (“Soil column 2” in Figure 2) since it most corresponds to the depth below the Station IV, and the simulated surface diffusive flux, which is by definition averaged over the reservoir surface (Figure 11).

The main task of using the model is to estimate the annual emission of methane, which implies a qualitative reproduction of seasonally averaged emission values, as well as significant methane emission due to the mixing of the water column during the destruction of stratification or hydrostatic pressure drops. The model satisfactorily reproduces both background values of the methane flux in the summer period and episodic emissions, which is clearly seen in 2018 and 2019 (Figure 11)—the extrema of the flux values according to observations correspond to the peaks according to the model; however, they are shifted in time by 1–2 days. The model results were less successful in the summer of 2017. The methane flux significantly increased during the summer, likely due to the accumulation of fresh organic matter at the top of sediments. At the same time, in the model, the temporal variation of methane production in sediments is caused only by temperature dynamics. In the future, the model will be improved by incorporating the variable content of labile organic matter in bottom sediments.

It is noteworthy that oxygen supersaturation in the epilimnion can be a “predictor” for a sharp summer increase in methane emissions since this supersaturation indicates significant primary production, which is not compensated for by respiration in the mixed layer, implying rapid removal of dead organic matter from the mixed layer, so that the matter, deposited at the bottom, becomes a substrate for methanogens.

The comparison between annual methane emission estimates by in situ measurements and the LAKE model is presented in

Table 8. Annual methane emissions from the Mozhaysk reservoir during 2015–2019 according to the LAKE model simulation and the estimation according to observations.

Year	Observations			LAKE Model
	Emission from Surface, tC-CH4	Degassing, tC-CH4	Total, tC-CH4	Diffusive, tC-CH4	Ebullition, tC-CH4	Degassing, tC-CH4	Total, tC-CH4
2015				50.0	287.5	6.5	344.0
2016				47.0	268.0	11.4	326.4
2017	334	13	347	51.7	273.8	11.8	337.2
2018	256	9	265	52.9	360.8	15.2	428.9
2019	378	11	389	60.5	354.9	8.9	424.3

.

The annual values of CH₄ fluxes to the atmosphere obtained in the basic run of the LAKE model and measurements are close, except for 2018. Most probably, this is related to temperature calculation in the LAKE model in the summer of 2018. Due to high bottom layer temperatures in the model, the methane generation rate exceeds the real methane generation rate. The average methane flux from the Mozhaysk reservoir to the atmosphere for the entire simulation period is 34.4 mg C-CH₄ m⁻² day⁻¹ according to field observations and 37.7 mg C-CH₄ m⁻² day⁻¹ according to the model. These values, when compared with global inventories, refer to the Mozhaysk reservoir as the upper limit of the emission range for reservoirs in the temperate zone [21].

The main contribution to the total methane emission according to the LAKE model results makes the ebullition flux, which is especially significant for a large total flux (as mentioned before—during the periods of maximum methane flux into the atmosphere, ebullition can reach up to 95% of contribution in total flux by model results) (Figure 10), as well as according to field observations. The tendency of total flux and bubble flux increase starts from the spring period to the end of summer and early autumn, with maximal emissions during the autumn convection phase clearly visible (Figure 10), which is also consistent with the field measurements. The increase in flux during the warm period in the model is primarily due to an increase in the temperature of bottom water layers and sediments, assuming that methane production is determined by the intensity of microbiological processes that directly depend on the temperature regime. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

3.3. Model Sensitivity Experiments

An additional model experiment was conducted to assess the sensitivity of the model to changes in the maximum methane oxidation rate in Michaelis–Menthen kinetics. The equation is following [61]:

$${O=V_{max}\frac{C_{CH4}}{C_{CH4}+k_{CH4}}}_{}$$

(11)

where O—methane oxidation rate, C_CH4—methane concentration, k_CH4—constant of half-saturation for methane oxidation, and V_max—potential maximum methane oxidation rate.

In this experiment, the maximal reaction rate was doubled from 0.1 mol day⁻¹ (as in baseline simulation) to 0.2 mol day⁻¹ (Figure 12). An increase in maximal oxidation rate led to a decrease in the mixed-layer CH₄ concentration, bringing it much closer to observed values.

Table 9. Difference in mean values and RMSE of methane surface concentration between model calculations and observations.

Year	Mean Delta, μmolCH4 L−1		RMSE, μmolCH4 L−1
	Vmax = 0.1	Vmax = 0.2	Vmax = 0.1	Vmax = 0.2
2016	0.48	0.30	0.50	0.32
2017	0.28	0.11	0.34	0.20
2018	0.13	−0.07	0.92	0.88
2019	0.16	−0.01	0.34	0.27
All years	0.21	0.02	0.64	0.59

shows the difference in mean values of methane surface concentration between model experiments and observations.

It is important to note that changes in methane oxidation equations as constant Vmax do not affect the ebullition flux because bubbles cannot become oxidized. A small part of rising bubbles in the water column can dissolve and then become oxidized, but this is an insignificant part of the flux.

The next experiment with model parameters was an experiment with the parameter q10 in the methane generation equation in bottom sediments. This parameter determines the methane generation multiplier—increasing methane emission for every 10 °C rise in temperature. The results of these experiments are shown in Figure 13. For this figure, methane bubble fluxes from Model Column 2 were taken (Figure 2) because the depths of this column were closer to the depth of Station IV, where measurements were carried out (Figure 1).

Parameter q10 is more significant for ebullition calculation because ebullition is the main contributor to methane emission into the atmosphere.

Table 10. Difference in mean values and RMSE of methane fluxes into the atmosphere between model calculations and observations.

Year	Mean Delta, mg C-CH4 m−2 day−1		RMSE, mg C-CH4 m−2 day−1
	q10 = 2.3	q10 = 3.0	q10 = 2.3	q10 = 3.0
2016	−1.8	32.3	53.76	58.49
2017	−57.7	−27.4	152.67	136.89
2018	49.6	98.2	177.06	214.52
2019	−125.8	−77.7	245.81	218.72
All years	−29.8	10.5	177.98	181.14

shows the difference in mean values between calculated and measured flux.

The results of annual methane emissions calculated from the Mozhaysk reservoir based on the results of these two experiments (V_max = 0.2; q10 = 3) in comparison with regular conditions (V_max = 0.1 and q10 = 2.3) are shown in

Table 11. Methane emission from the Mozhaysk reservoir based on the results of experiments with changes in the potential rate of methane oxidation and the methane generation parameter q10.

Year	Emission, tonC-CH4 Regular Conditions	Emission, tonC-CH4 Vmax = 0.2	Emission, tonC-CH4 q10 = 3
2015	344	318	443
2016	326	296	447
2017	337	306	430
2018	429	394	555
2019	424	393	528

.

These two parameters—V_max and q10 can be effectively used for model calibration: Vmax for calibration diffusive emission, q10—for ebullition. Calibration of these two parameters can give more accurate results of methane emission calculation.

4. Conclusions

A comparison of the LAKE2.3 model results with the observational data showed that the model reproduces the temperature regime of the Mozhaysk reservoir relatively well. The temperature of the bottom water horizons is reproduced significantly worse than the temperature of the epilimnion according to Pearson R and RMSE criteria, which is due to the complexity of the vertical heat transfer mechanisms under stable stratification in a reservoir with throughflow and suggests the development of new background diffusivity parameterization for one-dimensional models better representing the observed mixing. The dissolved oxygen content in the bottom horizons (a very important characteristic for simulation of methane accumulation in hypolimnion), according to the simulation results, is very close to the data obtained instrumentally.

The total annual methane emissions from the Mozhaysk reservoir based on field data and model calculations are close. The average flux during the simulation period was 34.4 mgC-CH₄ m⁻²day⁻¹ according to field observations and 37.7 mgC-CH₄ m⁻²day⁻¹ according to the model. This convergence of the results was achieved by calibrating the following parameters: the maximum potential rate of the methane oxidation (an effective parameter for regulating the diffusion flux) and the parameter for the temperature dependence of methane production q10 (shown to be effective for regulating the bubble flux). The bubble component makes the greatest contribution to the total methane flux. In the summer, it contributes 95% of the total flux. The contribution of diffusive, ebullition, and degassing fluxes averaged for 5 years of simulations are 52,309 and 11 tons of methane per year, respectively.

The LAKE2.3 model presented in this study can be calibrated using available measurement data at a given reservoir (exemplified in this study by Mozhaysk reservoir) and then used to reproduce the year-round variability of CH₄ flux including the periods where the measurements are not conducted (for example, during the methane release immediately after the ice-off) and thus providing better estimation of annual emission. The model can be used to simulate CH₄ emissions from reservoirs imposed on different water level management, as well as under contrasting scenarios of future climate change. Our numerical experiments demonstrated the key parameters of the model responsible for the realistic reproduction of methane fluxes at the water–atmosphere interface and better estimates of annual emissions from artificial reservoirs.