2.1 Modeling approach

In this study, we assessed the applicability of four different models for the gas transfer velocity, referred to as gas exchange models, to a process-based physico-biogeochemical lake model, MyLake C. The four gas exchange models were selected because their performance in estimating air–water CO₂ fluxes in a small boreal lake has been extensively assessed in previous studies by Heiskanen et al. (2014), Mammarella et al. (2015), and Erkkilä et al. (2018) by comparing the calculated fluxes with direct CO₂ flux measurements. The models include (1) the widely applied experimental wind-based regression formula by Cole and Caraco (1998), (2) a boundary-layer model developed by Heiskanen et al. (2014), (3) a surface renewal model by Tedford et al. (2014), and (4) a regression model by MacIntyre et al. (2010).

2.1.1 Parameterization of air–water gas exchange

The flux of CO₂ between water and the atmosphere, $F_{{\mathrm{CO}}_{\mathrm{2}}}$, can be parameterized as the product of the CO₂ concentration difference between the surface water and the atmosphere with the gas transfer velocity k (Cole and Caraco, 1998):

$$\begin{array}{}\text{(1)}& F_{{\mathrm{CO}}_{\mathrm{2}}}=\mathit{\alpha }k(C_{\mathrm{w}}-C_{\mathrm{eq}}),\end{array}$$

where C_w is the CO₂ concentration in the surface water below the air–water interface, C_eq is the equilibrium concentration of CO₂, that is, the water column CO₂ concentration in the state of equilibrium with the overlying atmosphere, and α is the chemical enhancement factor applicable for reactive gases, such as CO₂. Gas fluxes from water to the atmosphere are thus defined to be positive. If a lake is nonalkaline, α can be assumed to be 1 (Cole and Caraco, 1998). The equilibrium concentration is calculated by Henry's law as

$$\begin{array}{}\text{(2)}& C_{\mathrm{eq}}=K_{\mathrm{H}}\mathit{\chi }{p}_{\mathrm{a}},\end{array}$$

where K_H is the temperature-dependent aqueous-phase solubility (also known as the Henry's law constant) of CO₂ at surface water temperature, χ is the mole fraction of the gas in the atmosphere, and p_a is the atmospheric pressure.

The gas transfer velocity k can be simply parameterized as a function of wind speed alone, or more complex models can be applied to describe the air–water gas exchange process or the near-surface turbulence that governs the gas exchange. In each of the four gas exchange models assessed in this study, the parameterization of k is made using a different combination of parameters. The parameters of each model and their units are listed in Table 1. With the exception of the simple wind-based model by Cole and Caraco (1998), near-surface turbulence is driven in the models by both wind shear and thermal convection promoted by heat loss from the surface.

Cole and Caraco (1998)Heiskanen et al. (2014)MacIntyre et al. (2010)Tedford et al. (2014)

Table 1Parameters used in the parameterizations of the gas transfer velocity in the gas exchange models by Cole and Caraco (1998), Heiskanen et al. (2014), MacIntyre et al. (2010), and Tedford et al. (2014).

Download Print Version | Download XLSX

Convection-driven turbulence occurs when surface heat flux is directed out of the lake, that is, when the buoyancy flux is negative (MacIntyre et al., 2010). The buoyancy flux β is defined as (Imberger, 1985)

$$\begin{array}{}\text{(3)}& \mathit{\beta }={\displaystyle \frac{g{\mathit{\alpha }}_{\mathrm{w}}{Q}_{\mathrm{eff}}}{{\mathit{\rho }}_{\mathrm{w}}{c}_{p_{\mathrm{w}}}}},\end{array}$$

where g is the gravitational acceleration, α_w is the thermal expansion coefficient of water, Q_eff is the effective heat flux, ρ_w is the density of water, and $c_{p_{\mathrm{w}}}$ is the specific heat capacity of water. The effective heat flux is defined as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{Q}_{\mathrm{eff}}& =Q_{\mathrm{S}}+Q_{\mathrm{SW}}\left(\mathrm{0}\right)+Q_{\mathrm{SW}}\left(z_{\mathrm{AML}}\right)\\ & -{\displaystyle \frac{\mathrm{2}}{z_{\mathrm{AML}}}}\underset{\mathrm{0}}{\overset{z_{\mathrm{AML}}}{\int }}{Q}_{\mathrm{SW}}\left(z\right)\mathrm{d}z,\end{array}\end{array}$$

where $Q_{\mathrm{S}}=Q_{\mathrm{H}}+Q_{\mathrm{L}}+Q_{\mathrm{LW}}$ is the net surface heat flux, Q_H is sensible heat flux, Q_L is latent heat flux, Q_LW is net longwave radiation, Q_SW is shortwave radiation, and z_AML is the depth of the actively mixing layer (AML) (Imberger, 1985). All heat fluxes from the atmosphere into the lake are defined as positive. The last three terms in the equation represent the fraction of shortwave radiation that is trapped within the AML, denoted as Q_SW,AML. The attenuation of shortwave radiation at depth z in the water column can be calculated using the Beer–Lambert law:

$$\begin{array}{}\text{(5)}& Q_{\mathrm{SW}}\left(z\right)=Q_{\mathrm{SW}}\left(\mathrm{0}\right)e^{-K_{\mathrm{L}}z},\end{array}$$

where K_L is the total attenuation coefficient of shortwave radiation. The AML is defined as the near-surface layer in which the water column temperature is within a certain range, usually 0.02 ^∘C, of the temperature at the air–water interface (MacIntyre et al., 2001). The buoyancy flux is positive when the near-surface water is heating and negative under cooling conditions.

In the boundary-layer model developed by Heiskanen et al. (2014), near-surface turbulence is parameterized through wind-induced and convection-induced water-side velocity scales, which are characterized by the wind-induced water friction velocity at a reference depth, u_*ref, and the penetrative convection velocity w_*, respectively. The penetrative convection velocity is calculated as (Imberger, 1985)

$$\begin{array}{}\text{(6)}& w_{*}=(-\mathit{\beta }{z}_{\mathrm{AML}}{)}^{\mathrm{1}/\mathrm{3}}.\end{array}$$

The gas transfer velocity can also be parameterized by the total turbulent kinetic energy dissipation rate ε (MacIntyre et al., 1995). The rate can be measured directly or estimated from other measurable quantities with similarity scaling (Tedford et al., 2014). In the parameterization of ε by Tedford et al. (2014), both wind-induced stress and heat-induced convection generate turbulence near the lake surface during cooling, but wind is the only factor responsible for the turbulence during heating. The total turbulent kinetic energy dissipation rate is determined in terms of shear production ${\mathit{\epsilon }}_{\mathrm{s}}=u_{*\mathrm{w}}^{\mathrm{3}}/\mathit{\kappa }{z}^{\prime }$, where u_*w is the wind-induced water-side friction velocity, κ=0.4 is the von Kármán constant, and z^′ is a reference depth, and convective turbulence production ε_c, which equals the buoyancy flux β as

$$\begin{array}{}\text{(7)}& {\mathit{\epsilon }}_{\mathrm{TE}}=\left\{\begin{array}{ll}\mathrm{0.56}{\mathit{\epsilon }}_{\mathrm{s}}+\mathrm{0.77}\left|{\mathit{\epsilon }}_{\mathrm{c}}\right|& \text{if }\mathit{\beta }<\mathrm{0},\\ \mathrm{0.6}{\mathit{\epsilon }}_{\mathrm{s}}& \text{if }\mathit{\beta }\ge \mathrm{0},\end{array}\right.\end{array}$$

The wind-induced water friction velocity u_*w can be calculated from the atmospheric friction velocity $u_{*\mathrm{a}}=(\mathit{\tau }/{\mathit{\rho }}_{\mathrm{a}}{)}^{\mathrm{0.5}}$, where τ is the wind shear stress and ρ_a is the density of air, as in MacIntyre et al. (1995):

$$\begin{array}{}\text{(8)}& u_{*\mathrm{w}}=u_{*\mathrm{a}}({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}}{{\mathit{\rho }}_{\mathrm{w}}}}{)}^{\mathrm{0.5}}.\end{array}$$

2.2 Model application

We used the MyLake C application to Lake Kuivajärvi presented in Kiuru et al. (2018) as the basis of the study. The model setup, including model forcing data and the initial in-lake conditions, is nearly identical to that described in Kiuru et al. (2018). The minor differences are pointed out in Sect. 2.2.2.

2.2.4 Model calibration and validation

We estimated the MyLake C parameters utilizing a Markov chain Monte Carlo-based Bayesian inference algorithm following the procedures in the original calibration of the Lake Kuivajärvi application presented in Kiuru et al. (2018). Each of the four new versions of the MyLake C Lake Kuivajärvi application, using the models for k by Cole and Caraco (1998) (both the MyLake C version and the respective gas exchange model being hereinafter referred to as CC), Heiskanen et al. (2014) (HE), MacIntyre et al. (2010) (MI), and Tedford et al. (2014) (TE), was calibrated individually. The simulations with the MyLake C versions using different gas exchange models are hereinafter collectively referred to as GEMs. The model grid length was 0.5 m. The model was run from 8 January 2013 to 31 December 2014. The calibration period extended from 8 January to 31 December 2013, and the measurements in 2014 were used for model validation.

The calibrations were performed against the daily averages of the automatic water column CO₂ concentration measurements at the depths of 0.2, 2.5, and 7 m. We chose to apply the automatic measurements instead of the corresponding manual measurements used in the model calibration in Kiuru et al. (2018) because the calculation of daily CO₂ fluxes was based on the automatic measurements at 0.2 m in this study, and the simulation results were thus comparable with the calculated CO₂ fluxes. Even though the near-surface CO₂ concentration was the most significant factor considering air–water CO₂ exchange, deeper depths were included so that model behavior would also remain reasonable at deeper levels.

The calibrated model parameters were selected on the basis of the original calibration. However, because the new calibrations were not performed against water column DO concentrations, the parameters related to interactions between DO and CO₂, the photosynthetic quotient and the respiratory quotient, were excluded from the parameter set. The DIC inflow concentration scaling factor C_DI,IN, applied during open-water seasons, was introduced as a new calibration parameter. The other parameters included in the calibration were the vertical turbulent diffusion parameter a_k, the wind-sheltering coefficient W_str, the DOC-related specific attenuation coefficient of photosynthetically active radiation β_DOC, the maximal phytoplankton growth rate at 20 ^∘C ${\mathit{\mu }}_{\mathrm{20}}^{\prime }$, the phytoplankton death rate at 20 ^∘C m₂₀, the degradation rates of labile DOC k_DOC,1 and semilabile DOC k_DOC,2, the fragmentation rates of autochthonous POC k_POC,1 and allochthonous POC k_POC,2, and the sedimentary POC degradation rate k_POC,sed. The parameters obtained in the original calibration, or the default parameters, were used as the means of the prior parameter distributions.

One parameter chain with 3000 iterations was produced in each calibration. The starting points were set to 50th percentiles of the prior distributions. The first half of each resultant chain was discarded as a burn-in period, and the final parameters chains included 1500 parameter sets. The medians of the final posterior distributions (Figs. S1–S4 in the Supplement) were chosen as the calibrated parameters. They are presented, together with the default parameters, in Table 2. After the calibrations, additional goodness-of-fit metrics were calculated. The Nash–Sutcliffe (NS) efficiency gives a relative evaluation assessment, determining the relative magnitude of the residual variance compared to the variance of measurement data (Moriasi et al., 2007). The value of the normalized bias (B^*) describes a systematic overestimation ($B^{*}>\mathrm{0}$) or underestimation ($B^{*}<\mathrm{0}$) of a state variable in the simulation, whereas the normalized unbiased root mean square difference (${\mathrm{RMSD}}^{\prime *}$) shows if the standard deviation of the simulated values is higher (${\mathrm{RMSD}}^{\prime *}>\mathrm{0}$) or smaller (${\mathrm{RMSD}}^{\prime *}<\mathrm{0}$) than that of the measurements (Los and Blaas, 2010).

Table 2Calibrated model parameters for the different versions of MyLake C application to Lake Kuivajärvi with different incorporated gas exchange models (HE: Heiskanen et al., 2014, CC: Cole and Caraco, 1998, MI: MacIntyre et al., 2010, TE: Tedford et al., 2014). The default parameter values were used as the means of the prior parameter distributions.

Download Print Version | Download XLSX

3.1 Model calibration

Even though the differences between the formulations of the gas exchange models incorporated into MyLake C are rather notable, the resultant CO₂ concentrations did not differ substantially between the GEMs, that is, between the simulations with the MyLake C versions using different gas exchange models (Fig. 1). However, an optimal simulation result can be attained through many different combinations of processes related to in-lake carbon dynamics and fluvial and atmospheric exchange in MyLake C, which is seen in the variation between the parameter values obtained from the different calibrations (Table 2). The calibrations were performed only against CO₂ concentrations, and the aim of the calibration was not to try to reproduce the actual in-lake carbon cycling but rather to compare different possible ways to generate an optimal water column CO₂ concentration. The performance metrics for CO₂ concentration shown in the Supplement (Table S1) indicate that all GEMs yielded CO₂ concentrations ($B^{*}<\mathrm{0}$) that were too low at all depths during the calibration and validation periods with only a few exceptions. However, the CO₂ concentration measurements performed during the ice-covered periods were largely not applicable at 0.2 m because of the lake ice cover and sometimes also inapplicable at deeper levels because of incorrect functioning of the measurement system.

Figure 1Simulation results for CO₂ concentration with each GEM (mmol m⁻³) versus the daily averages of automatic high-frequency CO₂ concentration measurements at the depths of (a) 0.2 m, (b) 2.5 m, and (c) 7.0 m in Lake Kuivajärvi during the calibration year 2013 and the validation year 2014.

Download

The average near-surface (0–0.5 m) CO₂ concentrations over the open-water seasons were notably higher in CC (44.3 and 40.3 mmol m⁻³ in the calibration year 2013 and in the validation year 2014, respectively) than in the other GEMs (HE: 34.2 and 31.6 mmol m⁻³; MI: 31.5 and 29.4 mmol m⁻³; TE: 36.9 and 34.1 mmol m⁻³). Only the days with available corresponding water column CO₂ concentration measurement data were included in the averaging of the simulated near-surface CO₂ concentrations over the open-water seasons. By contrast, the averages of the measured near-surface (0.2 m) CO₂ concentrations over the open-water seasons were 45.2 mmol m⁻³ in 2013 and 37.2 mmol m⁻³ in 2014. Thus, CC yielded a higher near-surface CO₂ concentration compared to the measurements in 2014 when only the ice-free season, the period of air–water CO₂ exchange, is considered. The simulated open-water seasons were determined from the simulated ice-off and ice-on dates. Because CO₂ flux differs from zero starting from the day after ice-off in MyLake C, the simulated open-water seasons applied in the study were 3 May–25 November 2013 and 16 April–22 November 2014. In 2013, the observed open-water season lasted from 1 May to 27 November. In 2014, the observed ice-off date was 12 April.

The simulated CO₂ transfer velocities and air–water CO₂ fluxes are presented in Fig. S5. The yearly average values of k were lowest in CC and rather similar between the other GEMs (CC: 2.81 and 2.76 cm s⁻¹ for the calibration period and the validation period, respectively; HE: 5.44 and 5.33 cm s⁻¹; MI: 5.87 and 5.82 cm s⁻¹; TE: 4.73 and 4.66 cm s⁻¹). The differences in the simulated fluxes between GEMs were dissimilar to those in k because of the differences in the simulated near-surface CO₂ concentrations. The smallest k values in CC were compensated for by the highest near-surface CO₂ concentrations. By contrast, a high daily CO₂ efflux due to a high k in MI reduced the simulated near-surface CO₂ concentration compared to the other GEMs during the whole simulation period. Overall, the differences in yearly air–water CO₂ fluxes between GEMs were smaller than those in the values of k (CC: 0.22 and 0.20 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for the calibration period and the validation period, respectively; HE: 0.28 and 0.26 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$; MI: 0.25 and 0.24 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$; TE: 0.28 and 0.27 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$).

The CO₂ efflux during the first few days after ice-off was higher in GEMs with a high k, which increased the water column pH in comparison to CC. The differences remained rather constant during most of the open-water seasons. The near-surface pH was on average 0.20–0.26 and 0.18–0.25 units higher in the other GEMs than in CC during the open-water seasons of 2013 and 2014, respectively. As a result, the average fractions of CO₂ of DIC in the near-surface layer were respectively 6–8 and 5–6 percentage units higher in CC than in other GEMs, which also contributed to the higher near-surface CO₂ concentration in CC than in other GEMs. In addition, the open-water season average near-surface pH was 0.22 units higher in 2014 than in 2013 in all GEMs. Accumulation of bicarbonate in the water column in the course of the simulations may have resulted in an excessively high pH and thus a relatively lower CO₂ concentration in 2014 compared to 2013.

The differences in simulated temperatures between GEMs, primarily due to different attenuation of shortwave radiation in the water column, were rather small, especially at 0.2 m and at 2.5 m (Fig. S6). High epilimnetic concentrations of both Chl a and DOC, resulting from a low phytoplankton death rate and a high allochthonous POC fragmentation rate, respectively, in MI resulted in the strongest attenuation of shortwave radiation and thus the highest near-surface temperature because of a thinner and warmer epilimnion than in other GEMs. The open-water season average near-surface temperatures were 0.28–0.47 and 0.65–0.86 ^∘C lower than the corresponding measured averages in the calibration and validation periods, respectively, being highest in MI and lowest in TE. The differences were greatest in November before ice-on. The simulated near-surface temperatures tended to be somewhat too low in spring and early summer during both periods and somewhat too high in the late summer and autumn of the calibration year.

Heat transfer to the depth of 7 m right after the onset of the summer stratified period was insufficient in the calibration year in all GEMs, and small values of a_k also reduced heat transfer through the epilimnion during summer stratification. As a result, water column temperature remained too low at the depth of 7 m, which was located in the hypolimnion for most of the summer, during the stratified period in the simulations. However, the performance of the simulation of CO₂ concentration was also successful at the depth of 7 m. The summertime mixed layer thickness was rather similar between GEMs during the calibration year but more variable during the validation year. Simulated thermocline deepening matched the measurements during the late summer of the calibration year but was too early in the validation year. The deepening was slowest in HE because a somewhat stronger temperature gradient in the metalimnion, which was due to the smallest a_k, and resisted wind-induced thermocline erosion during summer.

3.2 Effective heat flux

The effective heat fluxes at the air–water interface simulated with each GEM on 3 May to 31 October 2013 and the corresponding values calculated on the basis of heat flux and radiation measurements are presented in Fig. 2a. The largest differences between the magnitudes and the directions of simulated and measured Q_eff were seen in early May. The simulated Q_eff was directed out of the lake throughout the study period except for a few occasions in early May and in October, whereas measurement-based calculations yielded more frequent occurrences of a positive daily Q_eff. Also, a negative Q_eff was often overestimated by the simulations because of overly high negative sensible and latent heat fluxes and net longwave radiation (Fig. S7). The performance of the simulation of the components of surface heat flux was rather poor (Table S2). Overall, the Q_eff simulation performance was not very good (R² = 0.39–0.41, RMSE = 48.2–49.2 W m⁻², NS = 0.11–0.14, $B^{*}=-\mathrm{0.47}\mathrm{\dots }-\mathrm{0.46}$, n=164). The differences in the simulated Q_eff between GEMs, resulting mainly from different surface temperatures, were quite small.

Figure 2(a) Daily effective surface heat fluxes (W m⁻²) simulated with each GEM and calculated on the basis of heat flux measurements. (b) Simulated and empirically determined depths of the daily actively mixing layer (m) in Lake Kuivajärvi in May–October 2013.

Download

The extent of shortwave radiative heating of the AML, Q_SW,AML, is dependent on z_AML. The simulated z_AML was greater than the measured daily average with a few exceptions at the beginning and near the end of the study period (Fig. 2b), which increased the simulated Q_SW,AML and decreased a negative Q_eff. The simulation with a daily time step generated clear temperature variation in the epilimnion only on days with a high amount of surface heating in early summer and midsummer, which resulted in an overly deep AML during most of the period. In addition, the model with a sequential description of thermal processes did not catch simultaneous wind mixing and surface heat exchange processes that resulted in a deeper observational AML in spring and late autumn. However, day-to-day variation in the discrepancy of Q_SW,AML was high throughout the study period. Also, the simulations highly underestimated the atmospheric friction velocity (R²=0.35, RMSE = 0.11 m s⁻¹, NS = −3.2, $B^{*}=-\mathrm{1.89}$, n=166) (Fig. S8), the simulated u_*a being on average only 46 % of the measured daily average. The simulated daily drag coefficient C_d at 1.5 m was affected by atmospheric stability conditions. The median C_d varied from $\mathrm{1.589}\times {\mathrm{10}}^{-\mathrm{3}}$ to $\mathrm{1.593}\times {\mathrm{10}}^{-\mathrm{3}}$ between the GEMs.

3.3 CO₂ exchange

The differences between simulated gas transfer velocities for CO₂ and the respective calculated values during the study period 3 May–31 October 2013 were rather small in the cases of gas exchange models based solely on wind speed, CC and MI, but the discrepancies were higher in HE and TE, which also include the effect of thermal convection on gas exchange (Fig. 3, Table S3). The simulations with CC and MI often yielded slightly higher values of k than the respective calculations because the simulated surface temperature was higher than the measured daily average (Fig. S6), and thus the temperature-dependent Schmidt number correction of k was different. Also, the occurrences of a simulated negative β in early May in MI yielded higher k_MI compared to the respective calculated values obtained from the observed positive β. The simulated k_HE was often higher than the calculated counterpart because of a high negative Q_eff or a deep AML in the simulations (Fig. 2), which resulted in a high penetrative convection velocity. In HE, the effects of wind-induced shear and thermal convection on k are set to be roughly of the same order of magnitude and the wind-induced shear velocity is calculated from wind speed, whereas CO₂ flux is driven principally by wind shear, which is calculated directly from u_*a, in TE. Because the simulated u_*a was consistently significantly lower than the corresponding daily measured average, the simulated k_TE was on average 40 % lower than the calculated value.

Figure 3Simulated and calculated gas transfer velocities for CO₂ (cm h⁻¹) in Lake Kuivajärvi on 3 May–31 October 2013 obtained with the gas exchange models by (a, e) Heiskanen et al. (2014), (b, f) Cole and Caraco (1998), (c, g) MacIntyre et al. (2010), and (d, h) Tedford et al. (2014).

Download

The simulated near-surface CO₂ concentrations were significantly too low during most of the study period in all GEMs except for CC, which yielded concentrations in autumn that were too high (Fig. 4, Table S3). The higher the simulated k and daily CO₂ efflux, the greater the resulting decrease in near-surface CO₂ concentration. The decline of the near-surface CO₂ concentration after ice-off was too rapid in all GEMs, especially in MI, the GEM with the highest k. The simulated near-surface CO₂ concentration also declined close to the atmospheric equilibrium concentration in all GEMs in late summer because of an insufficient gain of CO₂ in a shallow epilimnion developed under warm and calm conditions. The low simulated air–water CO₂ concentration gradients in May resulted in an underestimated air–water CO₂ flux from the water column compared to the respective calculated fluxes (Fig. 5, Table S3). The simulated flux was notably lower than the calculated flux in TE during the whole study period because of a small k_TE. In contrast, CC notably overestimated the corresponding calculated CO₂ flux in August and September because of a high simulated near-surface CO₂ concentration. Also, the simulated CO₂ flux was slightly higher than the calculated flux in HE in August and September because of high epilimnetic net CO₂ production. The total simulated CO₂ flux during May–October matched the calculated flux in CC but was notably lower in HE and MI and less than half of the calculated flux in TE (Table 3). The underestimated near-surface CO₂ concentrations were somewhat compensated for by the higher simulated k_HE and k_MI compared to the calculated counterparts, which decreased the difference between the simulated and calculated fluxes in HE and MI.

Figure 4Simulated CO₂ concentrations (mmol m⁻³) in the surface layer (0–0.5 m) obtained with each GEM and the daily averages of the automatic measurements at 0.2 m in Lake Kuivajärvi in May–October 2013. Also shown are the atmospheric equilibrium concentrations of CO₂ (C_eq) obtained from the simulations (dotted colored lines) and calculated from the measured atmospheric CO₂ concentration and surface water temperature (solid black line). Note the different vertical scales in May and in June–October.

Download

Figure 5Simulated and calculated air–water CO₂ fluxes ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) in Lake Kuivajärvi on 3 May–31 October 2013 obtained with the gas exchange models by (a, e) Heiskanen et al. (2014), (b, f) Cole and Caraco (1998), (c, g) MacIntyre et al. (2010), and (d, h) Tedford et al. (2014).

Download

Table 3Total and monthly averages of simulated and calculated CO₂ fluxes ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) in May–October 2013 obtained with different gas exchange models. Only the days with available measurement data are included in the averaging of the simulated fluxes. Monthly values for June are excluded because measurement data were available only for 7 d.

Download Print Version | Download XLSX

The applied gas exchange models yielded notably different calculated monthly CO₂ effluxes (Table 3). The CO₂ fluxes were calculated using the measured air–water CO₂ concentration gradients, and thus the differences between the calculated fluxes were only due to different values of k. Monthly fluxes calculated with MI were nearly or even more than double those calculated with the other wind-based model CC. Days with a positive β, resulting in a lower k_MI, occurred mainly in May and October, and thus the difference between the CO₂ fluxes calculated with MI and CC was slightly smaller in those months. The models that include the effect of thermal convection, HE and TE, yielded notably higher CO₂ fluxes than the simplest model, CC. Nevertheless, the CO₂ fluxes calculated with MI were slightly higher than those calculated with HE. The CO₂ fluxes calculated with TE were clearly the highest in all months, which was, however, not the case in the simulations.

The calculated daily values of k and CO₂ flux were dependent on the calculation interval. If the daily k had been calculated as the daily average of calculated half-hour values of k instead of using the daily averages of the input variables, the results would have been different. The daily averages of calculated half-hour k_MI (RMSE = 0.70 cm h⁻¹, $B^{*}=-\mathrm{0.16}$) and k_TE (RMSE = 0.22 cm h⁻¹, $B^{*}=-\mathrm{0.04}$) were lower than the respective values calculated using daily averages of input variables, whereas the opposite was the case for k_HE (RMSE = 0.48 cm h⁻¹, $B^{*}=\mathrm{0.20}$) and k_CC (RMSE = 0.16 cm h⁻¹, $B^{*}=\mathrm{0.15}$). In contrast, the calculation of a daily CO₂ flux as the average of half-hour fluxes yielded a slightly higher CO₂ flux in all GEMs (HE: RMSE = 0.066 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, $B^{*}=\mathrm{0.13}$; CC: RMSE = 0.034 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, $B^{*}=\mathrm{0.11}$; MI: RMSE = 0.10 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$,$B^{*}=\mathrm{3.4}\times {\mathrm{10}}^{-\mathrm{4}}$; TE: RMSE = 0.11 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, $B^{*}=\mathrm{0.05}$).

The differences resulting from the different methods of the calculation of a daily k can partly be explained by the behavior of the driving variables of the models. Using the daily averages of the input variables in the calculation may have smoothened out the effects of the spells of stronger negative buoyancy flux or a deeper AML that increase the half-hour k_HE and the effects of the occasions of positive buoyancy flux that decrease the half-hour k_MI. Daily averaging of wind speed may have cut out the rapid increase in k_CC under stronger wind conditions during the course of day due to the greater-than-linear dependence of k_CC on wind speed. By contrast, because the dependence of k_TE on u_*a is less than linear and the impact of thermal convection on k_TE is minor, the effect of the diel variation of u_*a and thus the relative difference between the methods of the calculation of k_TE was rather small.

3.4 Lake CO₂ budgets

The simulated CO₂ budgets for the epilimnion of the lake during periods of continuous summer stratification in 2013 and 2014 differed between GEMs as a response to different CO₂ effluxes (Table 4). The simulations were not able to reproduce the short-lived episodes of a very shallow epilimnion on days with high solar radiation and low wind speeds in late August and early September 2013, but at other times the simulated z_epi matched the depths estimated from the measured daily temperature profiles rather well (Fig. 6). The epilimnion formed 11 d earlier and extended to 7 m 16–22 d later in 2013 than in 2014. The in-lake CO₂ concentrations were higher at the onset of stratification in 2013 than in 2014 because of less effective water column ventilation during the shorter spring mixing period. As a result, the amount of CO₂ in the epilimnion decreased during the stratified period in 2013, whereas it increased slightly in 2014.

Figure 6Simulated and observed depths of the epilimnion (m) in Lake Kuivajärvi during the continuous summer stratification in 2013 and 2014. The simulations were performed using each of the gas exchange models incorporated into MyLake C.

Download

Table 4Simulated CO₂ budgets (kg CO₂) for the epilimnion of Lake Kuivajärvi during summer stratification in 2013 and 2014 using different gas exchange models incorporated into MyLake C.

^* The change in water column CO₂ content due to CO₂ efflux was approximately 1 % lower than the amount of CO₂ evaded because of consequent equilibrium reactions in the carbonate system.

Download Print Version | Download XLSX

A higher net in-lake CO₂ production or a higher terrestrial CO₂ load was required to compensate for the higher CO₂ efflux in GEMs that yielded higher values of k (Table 4). Phytoplankton concentration, regulated in MyLake C by growth and death rates μ^′ and m, respectively, impacts CO₂ dynamics both directly through the amount of carbon fixation and indirectly through changes in epilimnetic thermal structure due to the attenuation of solar radiation. A high μ^′ resulted in faster phytoplankton growth in spring and thus in an earlier occurrence of the spring bloom, and a small m resulted in a higher phytoplankton biomass during midsummer and late summer. HE and MI yielded the highest maximum near-surface Chl a concentrations, approximately 15 mg m⁻³ in 2013 and close to 20 mg m⁻³ in 2014. In CC and TE, Chl a concentrations were greater than 10 mg m⁻³ during the growth peaks but less than 5 mg m⁻³ at other times because of the high values of m. The open-water season average near-surface Chl a concentration was highest in HE (9.6 and 9.3 mg m⁻³ in 2013 and 2014, respectively), followed by MI (7.5 and 6.1 mg m⁻³), CC (3.9 and 4.0 mg m⁻³), and TE (2.3 and 2.1 mg m⁻³).

However, a high phytoplankton biomass did not imply high CO₂ consumption because of phosphorus limitation of phytoplankton growth in the model and the resultant reduction of photosynthetic CO₂ consumption under high Chl a and low bioavailable phosphorus concentrations in the simulations. Instead, CO₂ fixation occurred at a steady rate and the total CO₂ consumption over the whole growing season was relatively higher under a low Chl a concentration due to a high m. The highest average phytoplankton biomass in HE resulted in the highest CO₂ fixation; however, total net CO₂ production was also highest in HE because of high k_POC,1 and k_DOC,1. Small values of k_POC,1 and k_DOC,2 resulted in a relatively low net CO₂ production despite low CO₂ fixation and a high k_POC,sed in TE. Net CO₂ production was lowest in CC because of rather high total CO₂ fixation during the long growing season and the rather small k_POC,1 and k_DOC,1.

A considerable increase in the inflow DIC concentration by way of the scaling factor was essential in order to significantly increase the terrestrial CO₂ input to the lake in GEMs with a high CO₂ efflux. The measured inflow CO₂ concentration was 200–250 mmol m⁻³ until ice-off, less than 80 mmol m⁻³ during May, and mainly between 50 and 100 mmol m⁻³ during the summer and autumn. Thus, the default inflow CO₂ concentration was only approximately double the simulated near-surface CO₂ concentrations during most of the open-water season, and the effect of external CO₂ loading on in-lake CO₂ concentration was inevitably rather small, especially during the low-discharge period in late summer and autumn. The values of C_DI,IN determined the order of the amounts of the net external CO₂ load to the lake in the GEMs (TE: 42 000 and 45 000 kg CO₂ over the years 2013 and 2014, respectively; MI: 27 500 and 31 400 kg CO₂; HE: 26 500 and 30 600 kg CO₂; CC: 22 200 and 25 800 kg CO₂).

However, the total net external CO₂ loads to the lake over the stratification periods were slightly higher than the net external CO₂ loads to the epilimnion in Table 4 because stream inflow was directed into the metalimnion on days when the inflow temperature was lower than the epilimnetic temperature. The epilimnetic loads were 90 %–92 % and 98 %–99 % of the total loads in 2013 and 2014, respectively, the proportions being highest in CC and lowest in MI. The amount of CO₂ outflow was relatively large in CC because of the high epilimnetic CO₂ concentration; thus, the net external CO₂ load was relatively lower in CC than in other GEMs compared to the differences in C_DI,IN. In addition, because inflow pH was unaltered in the scaling of inflow DIC concentration, some of the increased CO₂ load was eventually evaded to the atmosphere in the simulations, but the bicarbonate fraction of DIC remained in the water column, which resulted in a slight increase in in-lake pH and a decline in the CO₂ fraction of DIC, especially in GEMs with a high k. Nevertheless, the impact of different amounts of bicarbonate loading on the in-lake pH was minor compared to the impact of different springtime CO₂ effluxes between GEMs.

4.1 Differences between calculated and simulated CO₂ fluxes

There was less variation between the air–water CO₂ fluxes simulated with different GEMs, that is, simulated with the MyLake C versions using different gas exchange models, than between the CO₂ fluxes calculated with the corresponding different gas exchange models on the basis on measured surface heat fluxes and air–water CO₂ concentration gradients (Table 3). This was caused both by differences between the simulated and calculated values of k and by insufficient epilimnetic CO₂ production in the simulations. An increased terrestrial CO₂ loading or an increased in-lake CO₂ production was needed to balance the higher CO₂ loss from the epilimnion through efflux in GEMs with a higher k compared to the simple wind-based CC (Table 4). Still, the simulations yielded near-surface CO₂ concentrations that were too low (Fig. 4, Table S3), which contributed to the underestimation of CO₂ fluxes (Fig. 5). Calibrating the model only against the near-surface CO₂ concentration and thus using even higher values for organic carbon fractionation and degradation parameters would have improved the performance of the simulation of epilimnetic CO₂ concentration; however, it would have resulted in uncontrollable and probably excessively high CO₂ concentrations in deeper layers, which is disadvantageous in a year-round, vertically layered lake model.

The day-to-day performance of the simulation of epilimnetic CO₂ concentration was also partly determined by the simulated thermal stratification and epilimnetic volume. The simulations generally yielded a near-surface CO₂ concentration that was too low when the simulated z_epi was in accordance with the observed depth and performed more adequately only during periods when the simulated z_epi was too high (Figs. 4 and 6). The measurements showed an increase in the near-surface CO₂ concentration when the epilimnion became thicker, and vice versa, during the stratified period in 2013. Thermocline tilting-induced upwelling and convection-induced entrainment transported more CO₂-rich water into the epilimnion on windy and cool days (Heiskanen et al., 2014). Conversely, high solar radiation input combined with calm conditions results in the warming of near-surface water and the formation of a thin epilimnion with a lower CO₂ concentration. High solar radiation also enhances photosynthesis and thus increases the uptake of CO₂ (Provenzale et al., 2018). An overly deep simulated epilimnion resulted in enhanced CO₂ release from deeper layers and a higher total net CO₂ production in a larger epilimnetic volume, which were able to compensate for the CO₂ efflux in the simulations.

The accuracy of the determination of a daily Q_eff and the applicability of the concept of a daily AML are issues that may cause uncertainties when gas exchange models are used either to calculate or to simulate daily estimates of k. The calculated half-hour Q_eff was generally directed into the lake on some occasions during daytime because of solar heating of the AML and always directed out of the lake at nighttime; z_AML often increased during nighttime and decreased under the radiative heating of near-surface water during daytime. Boundary-layer models and surface renewal models have been developed to describe the short-term dynamics of turbulence in a shallow AML, and thus they may not perform equally well in calculations with a daily time step.

The wind-based CC yielded the lowest and the surface renewal model TE the highest calculated air–water CO₂ fluxes, which is in line with the comparisons of different gas exchange models using data from Lake Kuivajärvi by Mammarella et al. (2015) and Erkkilä et al. (2018); however, the differences in simulated CO₂ fluxes between CC and other GEMs were notably smaller than the corresponding differences in the two experimental studies. The performance of TE is strongly dependent on the magnitude of u_*a because wind shear is highly dominant over thermal convection as the generator of turbulence in the model. Because the simulations yielded significantly lower u_*a compared to the values obtained through EC measurements (Fig. S8), the CO₂ flux obtained with TE was much lower than the corresponding calculated flux. Also, Erkkilä et al. (2018) found that u_*a calculated from wind speed was lower than the measured u_*a in Lake Kuivajärvi. Bulk models for surface stress may yield low values for u_*a over a lake, especially when parameterized for open-sea conditions with low surface roughness (Wang et al., 2015), which is the case in MyLake C. Lake size may also affect the relative differences between gas transfer velocities obtained with different gas exchange models. Dugan et al. (2016) applied different gas exchange models to the calculation of DO exchange in temperate lakes of various sizes. Simple, wind-based models yielded clearly lower values of k than more complex models in lakes similar to Lake Kuivajärvi in size, whereas the differences between the model types were smaller in larger lakes with generally higher wind speeds and a higher relative importance of wind-induced mixing compared to convection. In addition, ecosystem-specific empirical regression models may not be suitable for lakes with dissimilar characteristics (Vachon and Prairie, 2013).

4.2 Comparison to EC CO₂ flux measurements

Estimates of air–water CO₂ fluxes obtained with the gas exchange models applied in our study have been compared with 30 min block-averaged EC CO₂ flux measurements over Lake Kuivajärvi (Heiskanen et al., 2014; Mammarella et al., 2015; Erkkilä et al., 2018). Heiskanen et al. (2014) compared the half-hour k calculated with HE, CC, and MI with those obtained through EC measurements of CO₂ flux in August–November 2011. In the study, the average values of k_HE and k_MI were approximately 70 % of the corresponding measurement-based values, but the average k_CC was only about half of the average k_HE and k_MI. Erkkilä et al. (2018) compared the daily medians of EC CO₂ flux during a 2-week period in October 2014 with the daily median CO₂ fluxes calculated with CC, HE, and TE. The CO₂ fluxes obtained with HE and TE were 60 % of the EC CO₂ fluxes and approximately double the CO₂ fluxes obtained with CC. Overall, TE yielded the best correspondence with the EC fluxes. TE also outperformed CC in the comparison of half-hour CO₂ fluxes during the open-water periods of 2010 and 2011 in Mammarella et al. (2015). In our study, the best agreement with simulated and calculated CO₂ fluxes was found in CC, whereas TE yielded the lowest simulated fluxes in comparison to the corresponding calculated fluxes. Thus, none of the GEM outputs can be considered compatible with EC CO₂ fluxes, provided that the conclusions from the half-hour comparisons in the abovementioned studies can be extended to a daily scale.

The simulation results for the daily air–water CO₂ fluxes cannot be directly compared with EC data because the data coverage of EC flux measurements is often low. For example, the data coverages for CO₂ flux were 27 % and 37 % in Erkkilä et al. (2018) and Mammarella et al. (2015), respectively. Quality screening excludes much of the measurement data, and short-time system malfunction may cause significant data loss during long study periods. Daily average or median EC CO₂ flux may not be representative for the whole day because of the temporal bias of the measurements. EC flux measurements often tend to be inapplicable, especially at nighttime, because of flux nonstationarity during light winds and cooling (Heiskanen et al., 2014) or the advection of CO₂ from the surrounding forest (Erkkilä et al., 2018). EC CO₂ fluxes over boreal lakes are often enhanced at night by water-side convection (Podgrajsek et al., 2015) or because of a higher air–water CO₂ concentration gradient due to the absence of photosynthesis as a CO₂ sink (Erkkilä et al., 2018).

Both the calculated and the simulated values of k were determined by means of the platform data. They were thus suitable for comparison with each other but may not represent the average conditions over the lake and hence may not yield correct estimates of whole-lake CO₂ fluxes. Wind speed, u_*a, Q_H, and Q_L were measured at a single point on the platform, and the source area of the EC measurements of u_*a, Q_H, and Q_L ranges from 100 to 300 m along the wind direction over the lake (Mammarella et al., 2015). Thus, the values may not be representative for the whole lake. Wind speed and the resulting u_*a over lakes surrounded by forests are lower in sheltered nearshore areas than in the central zones of the lakes (Markfort et al., 2010). Sheltering affects the spatial variation of wind speed, especially in small lakes such as Lake Kuivajärvi. Because Q_H and Q_L are dependent on wind speed over the lake, they may also be higher at the center of the lake than in nearshore areas. Also, the estimation of u_*a, Q_H, and Q_L in the simulations was based on wind speed and other forcing data obtained from the single-point measurements, and the simulated values may have been overestimates of the spatial averages. However, despite the same measurement location, some disparities existed between the simulated and measured Q_H and Q_L. The differences may be in part attributed to an underestimation of surface heat fluxes by the EC method, which was seen, for example, in a study on energy balance over a small boreal lake by Nordbo et al. (2011) and also in Mammarella et al. (2015). The sum of the measured EC heat fluxes in Lake Kuivajärvi was on average 83 % and 79 % of available energy in 2010 and 2011, respectively, in Mammarella et al. (2015). In addition, the total relative random error of the EC measurements is generally around 10 % for both sensible heat flux and latent heat flux as estimated in Mammarella et al. (2015). Considerable spatial variability may also occur in near-surface water CO₂ concentration in small, shallow boreal lakes (Natchimuthu et al., 2017), which may result in further discrepancies in the estimates on whole-lake CO₂ flux obtained on the basis of gas exchange models or by using a vertical, horizontally integrated lake model.

4.3 Factors influencing the epilimnetic CO₂ budget

The model parameter sets obtained through calibration of the MyLake C applications using different incorporated gas exchange models were notably different from each other, thus emphasizing different processes related to carbon cycling within the water column or to carbon exchange with the surrounding terrestrial ecosystem or the atmosphere. However, considering the main objective of the study, which is the simulation of near-surface CO₂ concentration and air–water CO₂ flux, the different outcomes of the calibration processes, that is, the different model parameter sets, can be considered equally justified as they give insight on the diversity of biogeochemical processes that impact lacustrine CO₂ dynamics.

Phytoplankton is a significant factor in the lake CO₂ budget and the main driver of the diurnal variation of CO₂ concentration in Lake Kuivajärvi (Provenzale et al., 2018). In MyLake C, inorganic carbon is fixed by phytoplankton and carbon is stored in autochthonous organic matter within the water column or in bottom sediments until it is mineralized by bacteria. A relatively large portion of epilimnetic phytoplankton and dead autochthonous particulate organic matter sank from the epilimnion into deeper layers in MI because of the small values of m and k_POC,1. Production of CO₂ via the degradation of phytoplankton-originated organic matter, as well as the release of bioavailable phosphorus in the epilimnion through the mineralization of autochthonous organic matter, was also slow in MI because of a small k_DOC,1. As a result, the net production of CO₂ in the epilimnion was rather low in MI (Table 4) despite the relatively high simulated phytoplankton biomass. Overall, differences in total net CO₂ consumption by phytoplankton during the stratified period between GEMs were rather small despite the large variation in the simulated phytoplankton biomasses because of the phosphorus limitation of photosynthesis in GEMs with a high phytoplankton biomass and because of the variation in the length of the active growing season between GEMs.

The simulated Chl a concentrations were rather constant over the growing season with the exception of the substantial spring growth peaks in CC and TE. There are no data on Chl a concentration in Lake Kuivajärvi in 2013, but the Chl a concentration at 0–3 m was at its highest, 30–50 mg m⁻³, in mid-July and decreased to a level of less than 2 mg m⁻³ in late autumn in the years 2011–2012 (Heiskanen et al., 2015). The epilimnetic Chl a concentration is usually 3–5 mg m⁻³ during the growing season with diatom-induced peaks under cool conditions in spring and autumn (Provenzale et al., 2018). Thus, the GEMs with low near-surface Chl a concentrations, CC and TE, may have yielded better estimates of the overall phytoplankton biomass than HE and MI. However, the net consumption of CO₂ by phytoplankton was not only related to the amount of phytoplankton biomass. Nevertheless, none of the GEMs captured the supposed monthly variation of epilimnetic CO₂ concentration caused by the seasonal succession of phytoplankton.

A conspicuously high C_DI,IN was needed to balance the high CO₂ efflux in the GEMs with a high k (Table 2). The restriction of the scaling of the inflow DIC concentration to the open-water season was a rough way to increase the gain of epilimnetic CO₂, and the summertime inflow CO₂ concentrations may have been unnaturally high, especially in TE. However, the use of C_DI,IN can be thought as the inclusion of the input of CO₂ through groundwater seepage to the lake. In budget calculations, groundwater DIC load can be generally estimated by applying groundwater DIC flow as a percentage of stream DIC load (Chmiel et al., 2016). The amount of inflowing groundwater and its properties in Lake Kuivajärvi are unknown. However, in addition to inflow through minor inlet streams and surface runoff, especially during snowmelt in spring, groundwater seepage may contribute somewhat to the total lake inflow volume because the measured total outflow volume over the year 2013 was approximately double the inflowing volume via the main inlet stream. The CO₂ concentration in groundwater in southern Finland is around 700–900 mmol m⁻³ (Lahermo et al., 1990), which is about tenfold higher than the estimated average inflow CO₂ concentration in Lake Kuivajärvi over the stratified period in 2013, 86 mmol m⁻³, and well in line with the yearly average groundwater CO₂ concentration near a boreal stream determined by Leith et al. (2015). Thus, groundwater-derived CO₂ transport to the lake may also affect the water column CO₂ concentration.

The effect of CO₂ inputs through minor inlets or groundwater may be supported by the fact that the simulated near-surface CO₂ concentration decreased too fast in all GEMs after ice-off in May 2013, that is, during a period when the snowmelt-induced flow in minor inlet streams may be substantial and when the groundwater level is generally relatively high (Fig. 4). The simulated epilimnetic CO₂ sinks were rather small at that time because net CO₂ consumption by phytoplankton was low in cool water and because CO₂ efflux was relatively low because of a low air–water CO₂ concentration gradient. Labile, autochthonous DOC was absent in the epilimnion in the simulations, and the degradation of allochthonous DOC was slow under the relatively cold conditions in May. Despite the measured inflow CO₂ concentration being approximately twice the simulated epilimnetic CO₂ concentration and the scaled inflow CO₂ concentrations and terrestrial CO₂ loads being even higher, the decline of the epilimnetic CO₂ concentration was rapid in all GEMs. The high abundance of diatoms in Lake Kuivajärvi in spring may have resulted in a supply of easily degradable organic matter, but net primary production also consumed CO₂. Thus, substantial CO₂ loadings through surface runoff, minor inlet streams, or groundwater seepage could have been plausible additional sources of epilimnetic CO₂ in May, provided that the additional surface inflow was rich in CO₂. The impact of groundwater seepage is supported by a study on the carbon budget of a small boreal lake by Chmiel et al. (2016), in which a discrepancy between estimates of the gain and loss of inorganic carbon was explained by a possible underestimation of the impact of groundwater inflow.

4.4 Implications for lake modeling

None of the four MyLake C versions with different gas exchange models surpassed the other ones in the study because of the complex interplay between the near-surface water CO₂ concentration and air–water CO₂ flux in the simulations. A higher CO₂ efflux would have required a higher gain of CO₂ in the lake through in-lake CO₂ production or external loading of inorganic carbon, but the MyLake C versions with gas exchange models yielding a high k were not capable of increasing the CO₂ gain sufficiently. Hence, it is not a trivial task to judge which of the four gas exchange models is most suitable for integration into MyLake C or other coupled physical–biogeochemical lake models. However, several experimental studies (e.g., Jonsson et al., 2008; MacIntyre et al., 2010; Heiskanen et al., 2014) have shown that traditional, wind-based models often yield low CO₂ fluxes when compared to estimates based on direct measurements. Thus, it is recommended to strive to use the more sophisticated and probably more correct gas exchange models provided that the biogeochemical lake model can be made adaptable to higher CO₂ losses and that the parameters included in the more complex, turbulence-based models can be correctly simulated. This also means that further improvements related to the description of in-lake carbon processes in lake models and to the modeling or other means of estimation of external inorganic and organic carbon loading are still needed. Despite the challenges in using complex process-based models in the assessment of carbon cycling in lakes, modeling is an effective means to quantify underlying processes related to lacustrine CO₂ emissions and to study the development of lake ecosystems under changing conditions.