Supplement of Numerical study of the seasonal thermal and gas regimes of the largest artificial reservoir in western Europe using the LAKE 2.0 model

Abstract. The Alqueva reservoir (southeast of Portugal) is the largest artificial lake in western Europe and a strategic freshwater supply in the region. The reservoir is of scientific interest in terms of monitoring and maintaining the quality and quantity of water and its impact on the regional climate. To support these tasks, we conducted numerical studies of the thermal and gas regimes in the lake over the period from May 2017 to March 2019, supplemented by the data observed at the weather stations and floating platforms during the field campaign of the ALentejo Observation and Prediction (ALOP) system project. The 1D model, LAKE 2.0, was used for the numerical studies. Since it is highly versatile and can be adjusted to the specific features of the reservoir, this model is capable of simulating its thermodynamic and biogeochemical characteristics. Profiles and time series of water temperature, sensible and latent heat fluxes, and concentrations of CO2 and O2 reproduced by the LAKE 2.0 model were validated against the observed data and were compared to the thermodynamic simulation results obtained with the freshwater lake (FLake) model. The results demonstrated that both models captured the seasonal variations in water surface temperature and the internal thermal structure of the Alqueva reservoir well. The LAKE 2.0 model showed slightly better results and satisfactorily captured the seasonal gas regime.


1 Representation of biogeochemical processes in LAKE model 1.1 Governing equations for dissolved gases and organic carbon in a water column Evolution and vertical distribution of three dissolved gases are considered in the LAKE2.0 model, which are 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 contributes to concentration of dissolved inorganic carbon (DIC), C DIC = C CO 2 +C HCO − 3 +C CO 2− 3 , and it is the change of DIC that reflects the number of carbon atoms in CO 2 molecules added to (or lost by) a solution from (to) atmosphere, bubbles, respiring organisms or decaying organical matter (see Section 1.2).
In addition, the content of dissolved organic carbon (DOC), particulate organic carbon (both living, POCL, and dead, POCD) are calculated. POCL includes carbon atoms contained in phytoplankton and zooplankton.
The species listed above obey the following equation system: where ∂z are diffusion operators, w g is a sedimentation velocity of POCD particles. Equations (4)-(6) do not contain A, because they are not derived by horizontal averaging, but follow from assumption of horizontal homogeneity of respective biogeochemical variable. This is caused by uncertainty of estimating the flux of these substances at the sloping interface between water and sediments. The r.h.s of these equations represent diffusion (assuming k s = k s,t + k s,m with the same eddy diffusivity k s,t and molecular diffusivity k s,m for all species; molecular dissusivity is not included in POCL and POCD equations), sources and sinks due to the following processes: • dissolution/exsolution of gases at the bubble-water interface (B CH 4 , B O 2 and B CO2 ); • photosynthesis (P O 2 , P CO 2 , P P OCL ); • respiration (R O 2 , R CO 2 , R P OCL ); • biochemical oxygen demand in the water column (D O 2 , D CO 2 , D DOC , D P OCD ); • sedimentary oxygen demand (S O 2 , S CO 2 ); • methane aerobic oxidation in the water column • death of living species (D h,P OCL ) 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 (1)-(3) are expressed in mol/m 3 that allows for simple relations of sinks and sources in different equations based on stoichiometry of the respective reactions. Organic carbon variables DOC, POCL and POCD in (4)-(6) are molar concentrations of carbon atoms contained in these organic groups. Terms B CO 2 , P CO 2 , R CO 2 , D CO 2 , S CO 2 , O CO 2 in (3) possess "CO 2 " subscript because carbon atoms are supplied to or removed from DIC of a solution in a form of CO 2 .
In the following, the parameterizations of processes related to O 2 and CO 2 dynamics are described, whereas formulations for CH 4 processes are presented in (Stepanenko et al., 2016).
The formulations for photosynthesis, respiration, biochemical oxygen demand and sedimentary oxygen demand basically adopted from (Stefan and Fang, 1994) and (Hanson et al., 2004).

Carbonate equilibrium
Carbonate equilibrium means the equilibrium in the following reactions: Involving kinetic constants of these reactions yields, that the DIC reads Here, the constants are given by Arrhenius 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, where pH is an external parameter.
Carbon atoms are added or removed from carbonate equilibrium system in a form of CO 2 during respiration, photosynthesis and organic chemical and physical processes, hence the change of C DIC equals to number of CO 2 consumed or produced. This explains the sense of terms in equation (3). For obtaining CO 2 flux across bubble surface or CO 2 diffusive flux to the atmosphere, C CO 2 is needed and is calculated from (9).

Boundary conditions for dissolved gases in a water column
The top boundary condition (at the lake-atmosphere interface) for any dissolved gas concentration in the case of open water has the form: 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: 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: with the Schmidt number Sc(T ) having individual values for different gases and being temperature-dependent. The k 600 coefficient has been a subject of numerous studies, and a number concepts have been put forward to quantify it (Donelan and Wanninkhof, 2002). The proper computation of this coefficient should account for the effects of a number of factors such as turbulence in adjacent layers of water and air, wave development and breaking, cool skin dynamics. The surface renewal model (MacIntyre et al., 2010;Heiskanen et al., 2014), used in LAKE2.0 model, "integrates" those effects through the near-surface dissipation rate of turbulent kinetic energy: where ν w designates molecular viscosity of water, C 1,SR = 0.5 is an empirical constant. TKE dissipation rate is available directly from k − ǫ closure.
When a lake is covered by ice, F C = 0, which neglects contribution of diffusion through ice cracks.

Photosynthesis
The intensity of photosynthesis in terms of oxygen molecules production is expressed as: The denominator here serves co convert units in 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: so that C P is a value of P max at the reference temperature T = T 0 . The limitation of oxygen production by the available photosynthetically active radiation PAR (S P AR ) is given by the Haldane kinetics: The PAR intensity delivering maximum to a limiter L min (=1) is S P AR = C Lmin,1 C Lmin,2 . In the model, these coefficients are specified as (Stefan and Fang, 1994;Megard et al., 1984): with H(•) denoting a Heavyside function, and T 00 standing for another reference temperature. It is seen from (17), that L min → 0 if S P AR → 0 and S P AR → ∞, i.e. PAR ihnibits photosynthesis at both low and high values of its intensity. The PAR instensity S P AR is expressed in a number of photons per square meter per hour, so that: where H sec = 3600 s and S * P AR is PAR intensity in W/m 2 . The coefficient transforming from J to Einstein (Einstein is an energy of Avogadro number of photons), T J→Eins , is estimated assuming the uniform distribution of energy in PAR region, which yields: with N A , h P , c denoting Avogadro number, Planck constant and the light speed in vacuum, respectively, all in SI units.
The treatment of chlorophyll-a concentration ρ Chl−a is given in Section 1.10.
Finally, from the gross photosynthesis reaction: or, in a shortened form: we see that the carbon dioxide consumption equals oxygen production, i.e. P CO 2 = P O 2 .
Equation (22) also implies that P P OCL = P CO 2 .

Respiration
P. Hanson et al. (Hanson et al., 2004) assume, that respiration is performed by "living particles", i.e. POCL, only in epilimnion, and may be scaled by gross primary production (i.e., photosynthesis rate), R P OCL = α P OCL P P OCL , α P OCL = 0.8. In contrast, we assume that this process happens at all depths where enough oxygen in situ to be used in respiration is available, with the same scaling. Evidently,

Biochemical oxygen demand (BOD)
We treat biochemical oxygen demand as a consumption of oxygen during degradation of dead organic particles (POCD) D P OCD and dissolved organic carbon (DOC) D DOC , following (Hanson et al., 2004); they suggest that D P OCD = ρ P OCD /τ P OCD , D DOC = ρ DOC /τ DOC with time scales τ P OCD = 20D sec , τ DOC = 200D sec (D sec is a number of seconds in a day). Thus, the BOD rate is:

Sedimentary oxygen demand (SOD)
The sedimentary oxygen demand appears as a sink in (2) and in essence is the contribution of the vertical flux of O 2 at the lake's bottom to the horizontally averaged oxygen concentration: Basing on the argument that SOD is controlled by both diffusion (governed by Fickian law) and biochemical consumption (described by Michaelis-Menthen kinetics), (Walker and Snodgrass, 1986) derive: where µ β 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: The stoichiometry of SOD is assumed to be close to that of 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 (2).  Hanson et al. suggest exudation to be scaled with photosynthesis rate, E P OCL = β P OCL P P OCL , β P OCL = 0.03 and the death rate to be defined as D h,P OCL = ρ P OCL τ Dh , where time scale τ Dh ranges from 1.1D sec in hypolimnion to 33D sec in epilimnion.

Sedimentation of organic particles
In the current model version we use the Stokes sedimentation velocity below the mixed layer: and the high-Reynolds-number limit of this variable in the mixed layer. Here, ∆ = ρ p /ρ w0 − 1, ρ 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 (Song et al., 2008), and the density of organic particles as 1.25 g/cm 3 (Avnimelech et al., 2001).

Chlorophyll-a dynamics
The chlorophyll-a dynamics in the model follows a simple scheme suggested in (Stefan and Fang, 1994), where chlorophyll-a density is calculated as: where the active layer, H a , is a maximum value between mixed-layer depth, H M L , and the photic zone depth, H P Z . The mixed-layer depth is defined as the depth of maximum Brunt-Väisälä frequency, and the photic zone depth is estimated as the depth at which the PAR irradiance drops to 10% of its surface value. The mean chlorophyll-a concentration in the active layer, ρ 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, α, using Poole and Atkins formula (Poole and Atkins, 2009): where z SD is the Secchi disk depth and k P A = 1.7. The above chlorophyll-a scheme is identical to that of (Stefan and Fang, 1994), excepting for it does not take into account the annual cycle of ρ Chl−a,0 .  Figure S1: ML water temperature, original and modified model results and their errors when compared to the observations.