2.1 Study site and data

2.1.2 Observational data

The lake has been monitored by the Dorset Environmental Science Centre (DESC) for more than 40 years, thus providing a comprehensive database supporting the model experiments. The high-frequency and continuous observations used in this study have been collected during 5 years (14 July 2010–19 October 2015) with 10 min resolution on a monitoring buoy in the lake center (Fig. 1). The meteorological variables collected are air temperature, humidity at a height of 1.75 m above the water surface, wind speed at 1.75 m, and downwelling shortwave and longwave radiation. Pressure and precipitation are available from a land-based meteorological station 500 m away from the western lake shore. The high-frequency dataset further includes a vertical temperature profile (up to a depth of 27.1 m prior to the summer of 2012 and to a depth of 9.85 m thereafter, aggregated to a time step of 1 h), O₂ concentrations (at two depths: 1, 18 m, aggregated to a time step – 1 h) and CO₂ concentration (at a depth of 0.39 m, for the period 12 March 2012–19 October 2015, aggregated to a time step of 1 h). In addition, traditional temperature and O₂ profiles were collected to a depth of 35 m, except during the wintertime period, at a time step of approximately 2 weeks, and Secchi depth was also measured fortnightly or monthly during these 5 years (for more details, see Tables S6 and S7).

The meteorological conditions of the region during the period under study are typical for continental climate of mid-latitudes (Table 1). This region is subject to the influence of the Azores High during summer and is located in between the Icelandic Low and the Canadian High forming in the central part of the continent during winter. Cold arctic air usually invades this area during winter, leading to the formation of stable freezing weather. However, it should be considered that winters are usually characterized by several abrupt rises in air temperature above 0 ^∘C (see Fig. S1 in the Supplement). This phenomenon causes the 2012 winter to rank as the warmest winter during the 5 measurement years. Average annual wind speed is relatively low (≈2 m s⁻¹), likely caused by wind-shading effects imposed by the surrounding forest. Predominant wind directions are WSW (230–280^∘) and SE (130^∘).

Table 1Hydrometeorological conditions during the period 14 July 2010–19 October 2015.

^* Water temperature measured at 0.1 m depth.

Download Print Version | Download XLSX

2.2 Lake models

In this study, a suite of 1-D lake models was applied to Harp Lake. Two classes of finite difference models can be identified based on different approaches of calculating the vertical ED: Henderson-Sellers-based models (Henderson-Sellers, 1985; Hostetler and Bartlein, 1990), which are computationally inexpensive, and the more sophisticated models based on the k−ε turbulence closure scheme. The Henderson-Sellers approach computes the ED coefficient as a function of the Richardson number and considers wind-driven turbulence during stable or neutral stratification. A separate procedure is considered for convective mixing and assumes instantaneous mixing of unstably stratified depth layers. The k−ε formulation involves prognostic equations for turbulent kinetic energy and its dissipation. Both wind-induced and convective turbulence are taken into account in this approach. Within the first class there are the recently developed models ALBM and MTCR-1, while k−ε models include LAKE and LAKEoneD.

The FLake model stands aside from the other 1-D models due to the two-layer bulk structure which employs the concept of self-similarity to estimate the temperature profile in the mixed layer and thermocline, respectively (Mironov, 2008). In the mixed layer the temperature is assumed to be constant, whereas below it is parameterized as a function of non-dimensional depth. ED is not explicitly specified in this model. Of all models, FLake has the lowest computational cost (Thiery et al., 2014a).

Ice formation is one of the main controlling factors for the dissolved gas dynamics in a lake during winter. Models like LAKEoneD and MTCR-1 involve a simple ice parameterization, in which ice cover is governed only by air temperature and lake temperature is not considered (Ashton, 1980, 2011). In the FLake model, the temperature profile within the ice cover is parameterized via a time-varying shape function with a linear asymptotic. The ALBM incorporates a more comprehensive ice parameterization developed by Fang and Stefan (1994), taking into account the physical processes of freezing of the surface water, radiation and heat penetration through ice and snow. The LAKE model takes into account the same physical processes and uses a multilayer ice parameterization (Stepanenko and Lykossov, 2005; Stepanenko et al., 2011) based on an unsteady heat transfer equation. Light extinction in all models follows the Beer–Lambert law with a constant light attenuation coefficient.

The ALBM and LAKE incorporate the most comprehensive biogeochemical modules, including the sources, sinks and transport of CO₂, O₂, and CH₄ (and nitrogen in the ALBM) (see Table S1 in the Supplement). LAKEoneD model calculates O₂ concentration only. Parameters and characteristics of all models are summarized in Table 2. All the models use finite-difference numerical schemes along vertical coordinates, except for FLake, which is based on the two-layer representation of the lake temperature profile (see above). The vertical grid of these models is fixed in time, while the number of numerical layers varies from 51 to 105 between models. The time step varies within the range of 20–3600 s between the models. The ALBM uses a fourth-order adaptive Runge–Kutta–Fehlberg solver and a variable time step. If there are abrupt changes in boundary conditions (such as air temperature, solar radiation and wind), the time step in this model will be reduced to avoid numerical instability.

Table 2Characteristics of lake models relevant to the study.

Download Print Version | Download XLSX

2.3 Experimental setup

In this study, a set of numerical experiments was conducted to compare model simulation outputs assuming different morphometries, light extinction coefficients, and for the biogeochemical model parts initial CO₂ concentrations (Table 3):

1. a reference simulation in which the depth is set to a maximum value (37.5 m), light attenuation coefficient is set to a mean measured value (μ=0.4 m⁻¹) and the morphometry is not taken into account (dependence of horizontal cross-section area on depth (hypsometric curve) is not included) (RefSim);

2. model simulations testing the sensitivity of models to variations in

   • lake depth, setting it to mean value h_ave=13.32 m (LDSim),

   • the light attenuation coefficient, prescribing μ_min=0.28 m⁻¹ and μ_max=0.68 m⁻¹, which are minimal and maximal values measured, respectively (ExtMinSim, ExtMaxSim),

3. a model simulation testing the effect of the vertical diffusion on the distribution of the passive tracer (PassTrSim), and

4. a model simulation with biogeochemical modules activated including the calculation of the O₂ and CO₂ concentration (GasSim). Measured mean CO₂ concentration and O₂ observations are used as initial profiles. The time step of model output in all experiments is 1 h.

Table 3Parameters of numerical experiments of LakeMIP-Harp.

^* For more details, see the simulation protocol in SI. Some of the experiments are not included in our study.

Download Print Version | Download XLSX

The concept of the self-similarity in the FLake model (Sect. 2.2) includes the “shape factor” coefficient C_θ, which determines the vertical temperature profile below the mixed layer. Evolution of this parameter is controlled by a relaxation timescale t_rc (Mironov, 2008). This timescale includes the dimensionless relaxation constant C_rc with a default value of 0.003 in the model. This constant is a calibration parameter which is generally individual for each lake (Shatwell et al., 2016; Kirillin et al., 2017). We tested the sensitivity of the model to the variation of C_rc using values 0.03, 0.3, 0.5, 1, 2, 3, and 30 and achieved the best agreement with the measured temperature profile using C_rc=0.3 (for details, see Sect. 3.1). In the following sections, this C_rc setting is used in RefSim and other experiments with the FLake model.

The meteorological forcing is identical for all models and experiments. This, however, does not imply that heat fluxes to the atmosphere are the same in models, as models use different surface flux schemes. Some of the models do not include bathymetry (e.g., FLake); thus, for a more consistent comparison, in the RefSim experiment all models are integrated assuming homogeneous depth. In the GasSim experiment, lake morphometry is considered because the lake's carbon budget cannot be accurately calculated without including sediment oxygen demand along a sloping bottom.

The effect of the tributaries is neglected due to the large residence time of Harp Lake. Likewise, water-level variations are not considered in the models given the small temporal fluctuations.

The choice of the lake depth in the 1-D model without bathymetry may be crucial for the modeling results, not only for the temperature regime (Balsamo et al., 2010), but also for the gas dynamics. The gas dynamics is affected by lake depth as the latter controls the amount of bottom-originated gas to be biochemically transformed in the water column (e.g., oxidation of CH₄). Thus, it is important to examine the sensitivity of models to this parameter. In all experiments except LDSim, the maximal depth is used to compare the output with all available measurements. To test the effect of the chosen depth on model performance, the depth is set to mean depth (13.32 m) in the LDSim experiment; as such, the modeled lake volume is equal to the real lake volume.

The light attenuation coefficient is an important parameter, regulating the thermal regime of a lake (Heiskanen et al., 2015) and biogeochemical processes including photosynthesis. Secchi depth measurements during the 5-year study period on Harp Lake are irregular but enable calculation of mean, maximum and minimum light attenuation coefficients using the Poole and Atkins formula (Poole and Atkins, 1929).

The motivation for conducting experiments with perturbed input parameters (depth and light attenuation coefficient) and initial conditions (dissolved gases) is that these properties, especially the light attenuation coefficient, become considerably uncertain in global applications. While a first global lake depth database recently became available (Choulga et al., 2014; Kourzeneva, 2010), such a dataset currently does not exist for water transparency.

To understand the role of the ED coefficient in gas distribution and dynamics, an additional experiment solving the pure vertical diffusive transport equation for a passive tracer governed by eddy diffusivities resulting from different lake models was conducted. Details on this experiment are presented in Sect. 3.3.

Harp Lake is a deep lake suitable for studying the transport of O₂ and CO₂ as well as their biogeochemical transformations due to long-term regular measurements of these gases. In this respect, the additional numerical experiment (GasSim) was conducted with the ALBM and LAKE, the lake models with the most sophisticated biogeochemical modules. To harmonize the experimental setup, the contribution of CO₂ and O₂ fluxes through a sloping bottom to the lake carbon budget was switched on in the LAKE model. It should be noted that the ALBM includes a CO₂ sink/source due to outlets and inlets, whereas LAKE does not take this effect into account. In the ALBM, a set of biogeochemical parameters was calibrated, specifically, maximum chlorophyll-specific photosynthetic rate, maximum metabolic loss potential, aquatic dissolved organic matter (DOM) microbial degradation rate, terrestrial DOM microbial degradation rate, and groundwater dissolved organic carbon (DOC) concentration (Tan et al., 2017). No parameter calibration was performed for the LAKE model (a complete set of default values of model parameters is given in Stepanenko et al., 2016).

The model skill scores used in this study are listed in Table 4.

Figure 2Time–depth pattern of temperature in Harp Lake (14 July 2010–19 October 2015), RefSim and observed data. The grey boxes represent the duration of the ice-cover period, and the white patterns denote the absence of data.

Download

Table 4Statistical measures of model performance.

Download Print Version | Download XLSX

3.1 Temperature and ice

All the lake models simulate seasonal water temperature dynamics adequately (Fig. 2): the epilimnion is warming throughout the first half of summer, accompanied by thermocline deepening, eventually leading to overturn in autumn. During winter, the lake temperature is at the freezing point immediately below the ice cover, causing stable stratification. During spring, the lake starts to mix due to the penetration of shortwave radiation. The ice formation and ice melt dates play a key role in various processes: the interaction of the lake surface with the atmosphere, the overturning during spring, biogeochemical processes, and release of gases that have accumulated during wintertime, e.g., CO₂ or CH₄. Yet none of the models is able to precisely predict these dates. The difference between the model results and observations become even larger towards the end of the simulation, with differences exceeding 2 weeks in winter 2014–2015 (see Table S3).

3.1.1 FLake model

In RefSim we test the model with different values of C_rc (see Sect. 2.3 above). The best agreement with the observed temperature profile we get using C_rc= 0.3. In particular, the error (RMSE_c) for surface temperature reduces from 3.4 ^∘C (using a standard value of c_relax=0.003) to 2.5^∘ and becomes close to other models. At the depth of ∼6–10 m RMSE_c reduces from 5.7 ^∘C (using a standard value of c_relax=0.003) to 5.2^∘. In addition, RMSE_c reduces from 3.3–0.7 to 1.8–0.2 ^∘C at the depths 20–27 m.

FLake demonstrates the largest error among other models for the depth of the thermocline (Fig. 3a) (RMSE_c up to 5.2 ^∘C). In contrast, the model performance in simulation of the lake surface temperature is comparable with other models, and the error is only 2.5 ^∘C (see Fig. S2a), which points to the fact that it includes similar formulations for the surface heat balance to the other models.

Figure 3(a) Centered root-mean-square error of modeled temperature calculated for each individual depth; (b) the correlation coefficient between the simulated and observed temperatures at each depth.

Download

As highlighted in former studies (e.g., Stepanenko et al., 2010; Thiery et al., 2014b), this model overestimates the surface mixed-layer depth during summer, even under ice (Fig. 2b). The mixed-layer depth in the FLake model is 1.5 times larger than in observed data (5.6 m versus 3.6 m). Predictions of ice cover and ice off are characterized by large errors (see Table S3), with modeled dates more than 2 weeks off from the observations.

Remarkably, FLake demonstrates the largest sensitivity among all models to the lake depth – RMSE_c is reduced approximately 2-fold (from 2.5 to 1.28 ^∘C) (see Table S4, Figs. S3 and S4) when using a mean depth of 13.32 m instead of the maximal depth 37.5 m. The thickness of mixed-layer depth reduces up to 4.5 m and it has a better agreement with the observations in comparison with the RefSim experiment.

FLake showed a similar response in the sensitivity model test to variation of the light attenuation coefficient to other models. However, the difference in temperature is smaller than in other models: in both RefSim and ExtMinSim, ExtMaxSim is less than 4 ^∘C in comparison with 8 ^∘C (Fig. 4). It can be related to the weaker temperature gradient in the thermocline in comparison with other models (Fig. 2).

Figure 4Temperature difference between the ExtMaxSim and RefSim experiments in three lake models: (a) ALBM; (b) FLake; (c) LAKEoneD.

Download

3.1.4 LAKEoneD and LAKE models

Even though LAKE and LAKEoneD use the same turbulence closure, they represent the temperature dynamics in the lake with considerable differences (Fig. 2c and d). In particular, the mean depth of the mixed layer as simulated by LAKE is 1 m deeper compared to LAKEoneD results (4.8 and 3.9 m, respectively). LAKE demonstrates greater heating of the epilimnion (up to 30 ^∘C) than LAKEoneD (maximum value of 27.3 ^∘C). The homogeneous distribution of temperature during autumn mixing predicted by LAKEoneD reaches the depth of 15 m in most years, as opposed to LAKE (up to 35 m) and the observational data (at least up to 25 m, where the deepest sensor is located). Remarkably, only in these two models do periods of unstable ice cover occur. For example, the ice in LAKE is thin and often comes off during winter (more than five times in winter 2011–2012), leading to wind-driven surface layer mixing (see Table S3).

Figure 5Time–depth distribution of the ED in Harp Lake (14 July 2010–19 October 2015), RefSim model experiment.

Download

The sensitivity experiments with varying light attenuation coefficients (see Table S4) generally result in a similar response in all models. Increasing the light attenuation coefficient (i.e., reducing water transparency) leads to an upward shift of the thermocline position (Fig. 4); conversely, decreasing the attenuation coefficient causes a downward move of the thermocline. This leads to a temperature change at respective depths of up to 4 ^∘C in the FLake model and 8 ^∘C in the ALBM and LAKEoneD. Figure 4 contains only three models with different types of turbulent closure because the effect of varying the light attenuation coefficient is very similar in models of the same type (e.g., k−ϵ models LAKEoneD and LAKE).

It is important to note that all models have a maximum error, RMSE_c for temperature (compared to observations) in the thermocline (Fig. 3a), which agrees with a previous study (Stepanenko et al., 2014). Conversely, the correlation between simulated and observed temperature decreases with depth in all models (Fig. 3b), starting from 0.9–1.0 near the surface to 0.1–0.7 close to the bottom.

3.2 Eddy diffusivity

Eddy diffusivity (ED), K, is the variable controlling the vertical turbulent transport in lake models. This section examines the results of simulated ED values, thereby focusing on discrepancies between models caused by different turbulent closure schemes.

The simulated time–depth distribution of ED is shown in Fig. 5. EDs in ALBM and MTCR-1 do not reproduce spring and autumn whole-water-column mixing as the Henderson-Sellers diffusivity is not valid for unstable stratification. In contrast, k−ε models correctly predict the deepening of the mixed layer in the ED pattern during summertime towards a complete mixing during autumn.

Even models employing the same type of turbulence parameterization differ in simulated eddy diffusivities in the thermocline by orders of magnitude. This is caused by the different background values of ED added to those estimated by turbulence closures and mimicking mixing mechanisms not directly resolved in 1-D models such as internal wave breaking (e.g., Hondzo and Stefan, 1993). Apparently, the background diffusivity is switched off in the MTCR-1 model during ice cover.

The increase in ED near the bottom during the stratified period in the ALBM can be explained by the change from a stable stratification in the thermocline to a quasi-neutral stratification beneath.

A peculiar feature of LAKE is that it features high-intensity turbulence during much of the ice-cover period in the lower part of the water column that is apparently related to residual flow after the autumn overturn. In contrast, in LAKEoneD, turbulence dissipates shortly after the ice cover. Both of these models demonstrate a subsurface region of high ED values in winter which is not reproduced by other models. Production of turbulent kinetic energy (TKE) below the ice may be caused by both momentum transfer through partial ice cover (implemented in LAKE) and shortwave radiation penetration causing under-ice convection. The conspicuous difference between LAKE and LAKEoneD results is that springtime deep mixing does not develop in LAKEoneD. This may be linked to different constants and stability functions used in its k−ε closures or different buoyancy and momentum fluxes at the lake–atmosphere interface during this season.

It is worth noting that the significant difference between the modeled and observed dates of ice cover and ice off (see Table S3) leads to errors in the simulation of ED as water opens and interacts with the atmosphere.

The main features of ED distribution over seasons described in this section hold in experiment LDSim with a different lake depth value (13.32 m).

3.3 Effects of eddy diffusivity on the vertical distribution of a passive tracer

The vertical distribution of ED in lake models is one of the major drivers of vertical transport of dissolved gases and their emission to the atmosphere. To isolate the effect of the ED parameterization in the different models, we solve the diffusion equation for a passive tracer concentration C:

The emission rate of the passive tracer at the bottom is set constant:

while at the surface the flux is assumed to be large enough to make the surface concentration negligible compared to deep-water concentration, so that the boundary condition here becomes

Setting zero concentration at the surface is justified by the fact that surface concentration of gases of interest (CO₂, CH₄) is typically orders of magnitude less than bottom concentration. The total (molecular plus turbulent) diffusivity K is taken from different models, so that the solution C of Eqs. (1–3) is individual for each lake model.

Calculation of tracer diffusion was carried out using simulated values K(z, t) from the two Henderson-Sellers-based models – ALBM, MTCR-1, and two k−ε models – LAKE and LAKEoneD, covering the period of 5 years. In order to take into account vertical mixing under unstable stratification performed in ALBM and MTCR-1 and not expressed by Henderson-Sellers diffusivity, a value K=10⁵ m² s⁻¹ was applied for the convective layer in such cases.

The results of the numerical tracer experiment are presented in Fig. 6. All models demonstrate seasonal variability of the passive tracer concentration, including the near-bottom accumulation during stratified periods. Three models (ALBM, MTCR-1, and LAKE) perform the full overturn of the water column during autumn and spring, leading to the emission of almost the whole tracer storage to the atmosphere and forming a homogeneous concentration profile. In LAKEoneD, the vertical mixing does not reach the lake bottom during spring and autumn in certain years (e.g., 2011 and 2014). Hence, not all tracer amount is removed from the water column during overturn periods (see Fig. S5b) in this model. As a result, not only the concentration of the tracer near the bottom, but also the total tracer amount in the water column increases in LAKEoneD throughout the simulation period. However, this effect is a result of the assumption of a bathtub-shaped lake bathymetry and would be less when including the depth dependence of lake volume.

Figure 6Vertical distribution of passive tracer concentration in Harp Lake simulations, according to different models.

Download

The MTCR-1 model demonstrates short overturn periods similar to those of LAKEoneD, but the mixing reaches deeper. Hence, the passive tracer does not accumulate over time. Furthermore, ALBM exhibits a higher near-bottom concentration than MTCR-1 because the ED simulated by ALBM is 1 order of magnitude lower compared to MTCR-1 at 30–35 m depth (see Fig. 6).

To assess the effect of different eddy diffusivities on the surface tracer flux, we use the ratio of the surface flux, F_surf, to the bottom flux, F_bot, where the former is given as

with C_eq=0 the surface concentration, C_w the concentration at the first model output level below the surface (0.3 m), both in mol m⁻³, Δz=0.3 m, and K_surf the surface diffusivity coefficient.

The Henderson-Sellers-based models produce short periods with high fluxes, primarily during autumnal mixing (see Fig. S5a), and the drops of the integral concentration in the water column at the same time (see Fig. S5b). The flux spikes are several orders of magnitude higher compared to k−ε models. A reasonable explanation is that instantaneous convective mixing implemented in the ALBM and MTCR-1 models leads to almost the whole tracer content inside the convective layer being released to the atmosphere during a single convective event. All the models demonstrate much weaker tracer fluxes during spring overturn compared to those resulting from autumnal convection.

Lake models with k−ε closure produce substantial tracer fluxes during the stratified period in particular years (e.g., 2013 and 2014). The likely reason for this effect is rapid changes in mixed-layer depth, due to sudden wind forcing, entraining higher concentrations from the hypolimnion into the mixed layer.

The tracer diffusion experiment was also conducted at the reduced depth of 13.32 m (see Fig. S5c and d). Due to the lower depth, seasonal vertical mixing now extends over the whole water column in all models, providing more effective release of the substance. ALBM and MTCR-1 demonstrate episodic high fluxes during summer along with k−ε models. The spring mixing in LAKEoneD now reaches the bottom in some years; hence, the integral tracer amount evolves here in a similar way compared to all other models. The time-averaged tracer emission in LAKEoneD is thus the most sensitive to depth reduction. The mean surface flux in this model increases up to approximately 54 %, while for the ALBM, MTCR-1, and LAKE this increase is 4 %, 19 %, and 8 %, respectively. This can be primarily explained by the contribution of the increased fluxes during summer stratification and of the autumn mixing now reaching the bottom.

Figure 7Time–depth profiles of O₂ in Harp Lake (14 July 2010–19 October 2015) and model experiment including sediments and morphometry (GasSim) for (a) ALBM and (b) LAKE. White patterns in (c) indicate absence of data.

Download

Figure 8Time series of O₂ concentration at depths 1 and 18 m, observed and simulated.

Download

3.4 Harp Lake biogeochemistry

3.4.2 Carbon dioxide

The vertical distributions of CO₂ concentrations simulated by the ALBM and LAKE are shown in Fig. 9. The ALBM simulates higher concentrations of CO₂ compared to LAKE. This is likely because, in contrast to LAKE, the ALBM considers dissolved inorganic carbon (DIC) inflow from the catchment and explicitly treats DOC and particulate organic carbon (POC) dynamics (Tan et al., 2017). Parameter calibration applied in the ALBM brought CO₂ content closer to measurements in this model. Correlation with the observational data at a depth of 0.39 m is relatively small for both models (0.3), and RMSE_c is 4.5 and 2.7 mg L⁻¹ for LAKE and the ALBM, respectively.

Figure 9Time–depth profiles of CO₂ concentration in Harp Lake (17 August 2012–20 April 2015), experiment including modules of sediments and morphometry (GasSim).

Download

Both models reproduce the seasonality of the CO₂ concentrations: during summer CO₂ accumulates in the hypolimnion, while CO₂ is consumed in the mixed layer due to photosynthesis. For the same reason, the depth of the O₂ maximum (Fig. 7c) is consistent with the metalimnetic CO₂ minimum.

Vertical turbulent diffusion (Sect. 3.3) greatly affects CO₂ patterns in models. In particular, the mixing periods clearly seen in Fig. 9 coincide with those in Fig. 6. Accumulation of CO₂ below the mixed layer during stratified periods is clearly demonstrated by both models. However, the vertical CO₂ distribution in the thermocline is very different; i.e., in LAKE, concentration attains its maximum near the bottom, whereas in the ALBM the maximum is located at 15–20 m depth. This can be related to the fact that in the ALBM, DOC and DIC from surface water and groundwater can be injected in the middle of the water column.

In contrast to oxygen, measured CO₂ content is much less correlated with simulated series. The wintertime increase in surface concentration seen in calculated series is hardly observable in measured data. In summer, simulated CO₂ series are smooth, whereas there are significant fluctuations in empirical data. There is some correspondence between seasonal-mixing-caused peaks in modeled concentration and the measured one, though with time lags. All this leads to the conclusion that, compared to oxygen, carbon dioxide is more controlled by biogeochemical processes misrepresented in lake models than by physical factors. As stated in the previous section, a realistic decay of oxygen content during stratified periods in the ALBM and LAKE models suggests that CO₂ amount produced by aerobic decomposition of organic matter in both the water column and in the top part of sediments is reasonably simulated as well. Satisfactory agreement of computed oxygen in the mixed layer and below with observations implies that photosynthesis minus respiration rate is fairly captured as well, and so do the models for CO₂ gain or loss due to these processes. Hence, the primary drawback of the models used in respect to DIC simulation is likely to be not explicitly simulating transport of carbon species from catchment to a water body. Thus, modeling approaches coupling the catchment and a lake presented recently (Futter et al., 2008; Duffy et al., 2018; McCullough et al., 2018) should be elaborated and used more widely.

The important role of catchment processes in building up the DIC levels in lakes introduces an extra difficulty for implementing lacustrine CO₂ emissions in the Earth system models. This is caused by the necessity to provide regional or global data on a lake catchment's geometric, physical and biogeochemical properties, which are not currently available. In relation to this, we can also note a faster progress on a roadmap to introduce lake CH₄ dynamics in climate models, where simulations passed from site-level studies to regional estimates (Tan et al., 2015), whereas CO₂ modeling is currently confined to individual lakes (Stepanenko et al., 2016; Tan et al., 2017; Kiuru et al., 2018).