Applicability and consequences of the integration of alternative models for CO2 transfer velocity into a process-based lake model

Freshwater lakes are important in carbon cycling, especially in the boreal zone where many lakes are supersaturated with the greenhouse gas carbon dioxide (CO2) and emit it to the atmosphere, thus ventilating carbon originally fixed by the terrestrial system. The exchange of CO2 between water and the atmosphere is commonly estimated using simple wind-based parameterizations or models of gas transfer velocity (k). More complex surface renewal models, however, have been shown to yield more correct estimates of k in comparison with direct CO2 flux measurements. We incorporated four gas exchange models with different complexity into a vertical process-based physico-biochemical lake model, MyLake C, and assessed the performance and applicability of the alternative lake model versions to simulate air–water CO2 fluxes over a small boreal lake. None of the incorporated gas exchange models significantly outperformed the other models in the simulations in comparison to the measured near-surface CO2 concentrations or respective air– water CO2 fluxes calculated directly with the gas exchange models using measurement data as input. The use of more complex gas exchange models in the simulation, on the contrary, led to difficulties in obtaining a sufficient gain of CO2 in the water column and thus resulted in lower CO2 fluxes and water column CO2 concentrations compared to the respective measurement-based values. The inclusion of sophisticated and more correct models for air–water CO2 exchange in process-based lake models is crucial in efforts to properly assess lacustrine carbon budgets through model simulations in both single lakes and on a larger scale. However, finding higher estimates for both the internal and external sources of inorganic carbon in boreal lakes is important if improved knowledge of the magnitude of CO2 evasion from lakes is included in future studies on lake carbon budgets.


Introduction
The majority of inland waters, especially in the boreal zone, are supersaturated with carbon dioxide (CO 2 ) with concentrations that can exceed the equilibrium concentration by several times and are therefore net sources of carbon to the atmosphere (Cole et al., 1994;Algesten et al., 2014). The contribution of lakes to the global carbon budget is recognized to be substantial in comparison to the role of marine and terrestrial ecosystems as global carbon sinks, but quantitative estimates of the global 5 contribution of lakes and other inland waters show significant variation (Cole et al., 2007;Battin et al., 2009;Tranvik et al., 2009). Atmospheric CO 2 exchange between lakes and the atmosphere is one of the key processes needed to be determined in constructing carbon budgets of lakes and in evaluating the role of lakes in global carbon cycling.
The exchange of weakly soluble gases, like CO 2 and oxygen, across the air-water interface is often modeled as a boundarylayer process in which the gas flux is proportional to the gas concentration gradient at the interface. The proportionality factor 10 k is known as the gas transfer velocity. In many long-used models for the gas transfer velocity, or gas exchange models, k is parameterized as a function of wind speed alone (Wanninkhof 1992, Cole andCaraco 1998). However, direct measurements of air-water CO 2 exchange using the eddy covariance (EC) method (Jonsson et al., 2008;MacIntyre et al., 2010;Heiskanen et al., 2014) have resulted in higher estimates of k compared to wind-based gas exchange models. For weakly soluble gases, k depends mainly upon turbulence in near-surface water (Banerjee, 2007), which is not generated merely by wind. Near-surface turbulence 15 is initiated predominantly by wind shear and negative buoyancy flux related to thermal convection induced by surface heat loss (Imberger, 1985). Buoyancy flux is relatively more important in small, wind-sheltered lakes, and parameterizations of the gas transfer velocity that are based solely on wind speed may not be applicable under such conditions (Read et al., 2012).
Turbulence-driven gas exchange models have been shown to be well in accordance with in situ measurements of k (e.g., Zappa et al., 2007;Vachon et al., 2010). 20 In surface renewal models, k is calculated as a function of the turbulent kinetic energy dissipation rate ε, which provides an indication of the intensity of near-surface turbulence (MacIntyre et al., 1995). Kinetic energy dissipation can be due to viscous and thermal processes, and ε is thus dependent on wind shear and convective heat flux (Lombardo and Gregg, 1989). Wind shear is characterized by the wind-induced water-side friction velocity. The water-side friction velocity can be estimated from the atmospheric friction velocity, which can be measured directly  or calculated by bulk formulas 25 using meteorological variables (Fairall et al., 1996). Heat-induced turbulence is generated if the surface heat flux is directed out of the lake. If measurements of the components of surface heat flux are not available, they can also be estimated using bulk formulas (Fairall et al., 1996).
Global estimates of carbon emissions from lakes often use conservative estimates of CO 2 fluxes or models that yield potentially underestimated values for k leading to low estimates of CO 2 fluxes (e.g., Cole et al., 2007;Raymond et al., 2013). Thus, 30 revised estimates of lacustrine CO 2 emissions will require higher net ecosystem production in the land-based ecosystems of the terrestrial biosphere to close the global carbon balance (Battin et al., 2009). Many studies concerning modeling of lake carbon balance (e.g., Bade et al., 2004;McDonald et al., 2013) or determination of lake carbon budgets (e.g., Sobek et al., 2006;Stets et al., 2009;Chmiel et al., 2016) also use simple wind-based models for k. Potential subsequent underestimates in carbon efflux may have consequences for the interpretation of carbon budgets in single lakes . A higher efflux may result in a re-evaluation of the amount of net ecosystem production in lakes or it can mean that external carbon sources are inadequately accounted for in lake carbon budgets.
The efflux of CO 2 from a lake is sustained mainly by in-lake CO 2 production through bacterial or photochemical degradation of organic matter in water column or in sediment. Widely across the boreal zone, the importance of the degradation of 5 allochthonous organic matter as an inorganic carbon source in lakes is conspicuous (Jonsson et al., 2001;Sobek et al., 2003).
Also the direct loading of terrestrially produced dissolved inorganic carbon (DIC) through surface water and groundwater inflows may lead to high CO 2 concentrations in some lakes (Maberly et al., 2013;Weyhenmeyer et al., 2015;Einarsdóttir et al., 2017).
In this study, we evaluated the performance of different gas exchange models in the simulation of air-water CO 2 flux in 10 a boreal lake with a process-based lake model and the adaptability of the lake model application to different CO 2 losses via efflux. We also calculated CO 2 budgets for the epilimnion of the lake during summer stratification on the basis of the simulation results and assessed the relative importance of different biogeochemical processes on the epilimnetic CO 2 conditions. We incorporated four alternative gas exchange models into a vertical process-based physicobiogeochemical lake model for the simulation of year-round profiles of water temperature and CO 2 concentrations with a daily time step. We then applied the lake 15 model to a humic boreal lake located in southern Finland for the period 2013-2014, calibrating each of the resultant alternative lake model versions against high-frequency water column CO 2 concentration measurements. We compared the simulated gas transfer velocities and air-water CO 2 fluxes with those calculated with the gas exchange models on the basis of measurement data. The aims of our study are (i) to assess the applicability of gas exchange models of different complexity to a process-based lake model with a daily time step and (ii) to assess the implications of higher CO 2 efflux estimates for the lake carbon budget. 20 2 Materials and methods

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 physicobiogeochemical lake model MyLake C. The four gas exchange models were selected because their performance in estimating air-water CO 2 fluxes in a small boreal lake has been extensively assessed in previous 25 studies by Heiskanen et al. (2014), Mammarella et al. (2015), and Erkkilä et al. (2018) by comparing the calculated fluxes with direct CO 2 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).

Parameterization of air-water gas exchange
The flux of CO 2 between water and the atmosphere, F CO2 , can be parameterized as the product of the CO 2 concentration difference between the surface water and the atmosphere and the gas transfer velocity k (Cole and Caraco, 1998): where C w is the CO 2 concentration in the surface water below the air-water interface, C eq is the equilibrium concentration 5 of CO 2 , that is, the water column CO 2 concentration in the state of equilibrium with the overlying atmosphere, and α is the chemical enhancement factor applicable for reactive gases, such as CO 2 . 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 10 where K H is the temperature-dependent aqueous-phase solubility (also known as the Henry's law constant) of CO 2 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 15 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.
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) 20 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 pw is the specific heat capacity of water. The effective heat flux is defined as where Q S = Q H + Q L + Q LW is the net surface heat flux, Q H is sensible heat flux, Q L is latent heat flux, Q LW is net longwave 25 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 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 30 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 windinduced and convection-induced water-side velocity scales, which are characterized by the wind-induced water friction velocity 5 at a reference depth, u * ref , and the penetrative convection velocity w * , respectively. The penetrative convection velocity is calculated as (Imberger, 1985) w * = (−βz AML ) 1/3 .
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.,10 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 ε s = u 3 * w /κz , where u * w is the windinduced water-side friction velocity, κ = 0.4 is the von Kármán constant, and z is a reference depth, and convective turbulence production ε c equaling the buoyancy flux β, as The wind-induced water friction velocity u * w can be calculated from the atmospheric friction velocity u * a = (τ /ρ a ) 0.5 , where τ is the wind shear stress and ρ a is the density of air, as in MacIntyre et al. (1995) 2.1.2 Gas exchange models 20 The widely applied experimental wind-based regression formula for k by Cole and Caraco (1998) gives the gas transfer velocity in units of cm h −1 as k CC = (2.07 + 0.215U 1.7 10 ) Sc 600 where U 10 (m s −1 ) is the wind speed at 10 m and Sc is the temperature-dependent Schmidt number of CO 2 .
In the boundary layer model by Heiskanen et al. (2014), the wind-induced water friction velocity is approximated to be a 25 linear function of the wind speed at 1.5 m height, U 1.5 : where C 1 is an empirical dimensionless constant, and the equation for k HE (m s −1 ) is where C 1 = 1.5×10 −4 and C 2 = 0.07 is another experimental dimensionless constant. The model by Heiskanen et al. (2014) is used in the vertical process-based Arctic Lake Biogeochemistry Model (ALBM) (Tan et al., 2017), which simulates inorganic and organic carbon cycling in permafrost lakes. The model by Heiskanen et al. (2014) is also included in the LakeMetabolizer package (Winslow et al., 2016), in which several lake metabolism models can be combined with models for computing the gas transfer velocity. 5 In the simple wind-based regression model by MacIntyre et al. (2010), the gas transfer velocity k MI (cm h −1 ) is calculated separately for heating and cooling conditions as In the surface renewal model of air-water gas exchange, k is parameterized as a function of the total turbulent kinetic energy dissipation rate as k = c(νε) 0.25 Sc −0.5 , where c is an empirical dimensionless constant and ν is the kinematic viscosity 10 of water (MacIntyre et al., 1995). Tedford et al. (2014) integrated the parameterization of the total turbulent kinetic energy dissipation rate, ε TE , into the surface renewal model to yield a model for the gas transfer velocity in units of m s −1 : The models by Cole and Caraco (1998)

Lake model MyLake C
We used an application of a one-dimensional process-based lake model MyLake C (Kiuru et al., 2018) for the simulation of the vertical distributions of water column temperature and CO 2 concentration and air-water CO 2 flux in the study lake. In addition, we integrated three alternative models for the gas transfer velocity into the lake model. MyLake C simulates inorganic and organic carbon cycling in a lake, taking into account terrestrial carbon loading, air-water exchange of CO 2 , and changes 20 in water column pH. However, groundwater exchange and changes in water level due to rainfall or evaporation are excluded.
The model operates on a daily time step, and the vertical grid length can be defined by the user. The model is based on a lake model MyLake v.1.2 (Saloranta and Andersen, 2007), which simulates lake thermal structure, seasonal ice and snow cover, and phosphorus-phytoplankton dynamics. In the model, vertical heat and mass diffusion are calculated with a diffusion equation using a vertical turbulent diffusion coefficient derived from the buoyancy frequency and parameterized by lake surface area by 25 default. Settling of particulate substances is also taken into account in the equation. In addition, convective and wind-induced water column mixing processes are included. As an exception to the daily time step, heat exchange between the water column and the atmosphere is calculated separately for daytime and nighttime. MyLake v.1.2 and its various extensions have been used in studies on stratification and lake ice cover (e.g., Saloranta et al., 2009;Dibike et al., 2012;Gebre et al., 2014), total phosphorus concentration and phytoplankton biomass (e.g., Romarheim et al., 2015;Couture et al., 2018), dissolved organic 30 carbon (DOC) concentration (Holmberg et al., 2014;de Wit et al., 2018), and dissolved oxygen (DO) conditions (Couture et al., 2015).
MyLake C has been designed to include only the most substantial physical, chemical, and biological processes related to carbon cycling in a well-balanced and robust way. CO 2 is produced in the lake through organic carbon degradation both within the water column and in the sediment and through phytoplankton respiration. Inorganic carbon production is coupled 5 to DO consumption, and vice versa. A division is made between readily degradable, phytoplankton-originated autochthonous particulate organic carbon (POC) and more refractory allochthonous POC. The model includes also the sedimentation, the resuspension and the permanent burial of POC. Correspondingly, DOC is classified into three compound classes with different bacterial degradabilities. A separate submodule (Holmberg et al., 2014) calculates the conversion of DOC into DIC via bacterial and photochemical degradation. The meteorological model forcing includes daily global radiation, cloud cover fraction, 10 atmospheric temperature, relative humidity, atmospheric pressure, wind speed at 10 m height, and precipitation. Hydrological forcing data include daily inflow volumes, inflow temperatures, inflow pH, and the inflow concentrations of modeled substances, including DOC, POC, and DIC. Complete data requirements are presented and model structure and applied equations are described in detail in Kiuru et al. (2018).
MyLake uses the Air-Sea Toolbox (Air-Sea, 1999) based on the parameterizations and algorithms in Fairall et al. (1996) for 15 calculation of surface wind stress and the components of surface heat flux. The sensible heat flux Q H , the latent heat flux Q L , and the wind shear stress τ are obtained from aerodynamic bulk formulas of the form 20 where c pa is the specific heat capacity of air, C h and C l are the transfer coefficients of sensible and latent heat, respectively, C d is the drag coefficient, U is wind speed, T a is air temperature, T s is water surface temperature, L e is the latent heat of evaporation of water, q a is the specific humidity, and q s is the saturation specific humidity at the water surface temperature. No wind sheltering effect on U is applied in the calculation of surface wind stress and surface heat flux components.
The air-water CO 2 flux F CO2 , given in units of mg m −2 d −1 in MyLake C, is calculated with Eq. (1) using the model for 25 k by Cole and Caraco (1998) (Eq. (9)). The chemical enhancement factor α is set to 1, and the temperature dependence of the aqueous-phase solubility K H is calculated according to Weiss (1974).
In this study, we incorporated the models for k by  (13)) into MyLake C as alternatives to the default model by Cole and Caraco (1998). The constants in the model by Tedford et al. (2014) are defined as c = 0.5 and z = 0.15 m as in Erkkilä et al. (2018). In MyLake C, the actively 30 mixing layer includes the model grid layers in which the water column temperature is within 0.02 • C of the temperature of the topmost grid layer. The temperature dependence of Sc for CO 2 is determined for surface water conditions using the polynomial fit in Wanninkhof (1992). The approximation U 10 /U 1.5 = 1.22 is used for the wind speed at different heights.

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).  (2005)). The length of the lake is 2.6 km, the maximum width is 0.3 km, and the surface area is 0.63 km 2 . The north-south-oriented lake has two distinct basins. The maximum depth of the deeper southern basin is 13.2 m (Heiskanen et al., 2014), which is more than double the mean depth 6.3 m. A measurement platform (Lake-SMEAR) is situated close to the deepest region of the lake. 10 The approximate retention time of the lake is 0.65 years. Lake Kuivajärvi is surrounded by managed mixed coniferous forest together with small open wetland areas . The majority of the catchment area (9.4 km 2 ) of the lake is flat. The main inlet stream with a mean pH of 6.5 (Dinsmore et al., 2013) drains four upstream lakes, which are smaller in area than Lake Kuivajärvi. The lake is dimictic: the spring turnover usually occurs rapidly right after ice-off in late April or early May, and the summer stratification period lasts until the autumn turnover in September or October. The duration of 15 the ice-covered period and the concomitant inverse stratification is usually 5-6 months . The turnover periods are hot moments for the release of CO 2 accumulated in the hypolimnion of the lake during stratification . Because of high terrestrial inputs of organic matter, a median concentration of DOC in the surface water is 12-14 mg L −1  and water clarity is rather low, a median light attenuation coefficient K L being around 0.6 m −1 . 20

Model forcing and calibration data
The meteorological forcing data and hydrological loading data used in the model application are described in detail in Kiuru et al. (2018). The daily averages of wind speed at 1.5 m and incoming shortwave radiation together with in-lake temperature and CO 2 concentration were obtained from automatic platform measurements (Heiskanen et al., 2014;Mammarella et al., 2015), and the remaining meteorological forcing data were obtained from the SMEAR II station or from weather stations 25 (Finnish Meteorological Institute) in Hyytiälä located less than 1 km from the lake (precipitation) and in Tikkakoski located approximately 95 km to the north-east from the lake (cloud cover fraction). Differing from Kiuru et al. (2018), the CO 2 mixing ratio in the atmosphere was assumed to be 395 ppm on the basis of the rather fragmentary time series of the high-frequency in situ measurements of the CO 2 mixing ratio, the method of which is described in Erkkilä et al. (2018).
The construction of the time series for lake inflow was based on continuous measurements of the discharges at the main inlet 30 and at the outlet of Lake Kuivajärvi in 2013-2014 (Dinsmore et al., 2013). Because the total measured outflow volumes were approximately double the main inlet discharge volumes on an annual scale, the daily inflow volumes were corrected by a factor of 2 in order to include the potential contributions of smaller inlet streams and groundwater to lake inflow. At the main inlet, water temperature was measured approximately two times a month in 2013 and continuously in 2014 and CO 2 concentration was measured two times a month in 2013 but mostly at intervals of 2-3 d around the period of ice-off in April and May using the procedure described in Miettinen et al. (2015). Daily time series were generated by linear interpolation.
The model was calibrated against the daily averages of the automatic high-frequency CO 2 concentration measurements: 5 an optimal set of selected model parameters were estimated so that the simulated CO 2 concentration time series matched the corresponding measured CO 2 concentration time series as well as possible. The estimation was performed using a statistical inference algorithm. In addition, the automatic water column temperature measurements were used in model performance validation. The CO 2 concentrations were measured at 0.2, 1.5, 2.5, and 7.0 m, and the temperature measurements were performed at 0.2 m, at 0.5 m intervals from 0.5 to 5.0 m, and at 6, 7, 8, 10 and 12 m using the measurement systems described in 10 Heiskanen et al. (2014) and Mammarella et al. (2015).

Model assessment data
We used additional meteorological measurements in assessing the performance of the alternative models for k incorporated into MyLake C during the period May-October 2013. An EC system located on the measurement platform measures the turbulent fluxes of momentum, heat, and water vapor (H 2 O) over the lake . The EC flux measurement system 15 includes an ultrasonic anemometer (USA-1, Metek GmbH, Elmshorn, Germany) and a closed-path infrared gas analyzer (LI-7200, LICOR Inc., Nebraska, USA) for measuring CO 2 and H 2 O mixing ratios at 1.8 m height above the lake surface. Air calculated from the H 2 O flux, and the atmospheric friction velocity was derived from the momentum flux. All EC measurement data were given as half-hour block averages. The EC measurements are explained in more detail in Erkkilä et al. (2018), and the description of EC data post-processing is found in Mammarella et al. (2015) and Mammarella et al. (2016). Contrary to the model forcing data, the air temperatures that were used in the measurement-based determination of the gas transfer velocities 25 were obtained from the platform measurements instead of the SMEAR II station when platform measurements were available.
In addition, the rather intermittent platform measurement data on relative humidity were used. In the calculation of the waterside friction velocity, missing relative humidities were replaced by a value of 75 %, which is close to the average of the SMEAR II measurements of relative humidity in May-October 2013, 72 %. The corresponding averages over the period May-August 2013, for which platform measurements were rather well applicable, were 66 % and 68 % for the SMEAR II and platform 30 measurements, respectively. Thus, the relative humidity can be assumed to have been slightly higher over the lake than at the SMEAR II station.
The estimation of the flux footprint distribution functions was made using the model by Kormann and Meixner (2001). The average footprint contributing to 80 % of the fluxes varies from 100 m up to about 300 m from the measurement platform depending on atmospheric stability conditions as described in Mammarella et al. (2015). Only wind directions along the lake (130 • -180 • and 320 • -350 • ) were included in the calculations to ensure that heat fluxes from the surrounding land were excluded. Furthermore, possible remaining effects of transversal advection during calm nights were removed through EC quality screening. In addition to the exclusion of some of the EC measurement data through the application of the quality screening criteria presented in Erkkilä et al. (2018), there was a gap in the heat flux data on 14-27 June because of EC 5 system malfunction. The monthly data coverage was 43-69 % and 32-70 % of the original data for sensible and latent heat fluxes, respectively. We constructed gap-filled half-hour time series for sensible and latent heat fluxes using linear fits between the measured sensible heat flux and wind speed multiplied by the air-surface water temperature difference and between the measured latent heat flux and wind speed multiplied by the vapor pressure difference, according to Mammarella et al. (2015).
Only the vapor pressures calculated from the measured relative humidities were used in the latter fit. The fitting was performed 10 independently for each month.
We compared the simulated gas transfer velocities for CO 2 and the simulated air-water CO 2 fluxes to those determined directly from measurements using the corresponding gas exchange models. The latter are hereinafter referred as to calculated gas transfer velocities and calculated CO 2 fluxes. The calculated CO 2 transfer velocities for each of the four gas exchange models were obtained using the daily averages of required measured variables. The calculated air-water CO 2 fluxes were 15 further obtained as the product of the calculated CO 2 transfer velocities and the daily averages of the measured air-water CO 2 concentration gradient. The conditions were thus compatible with the daily time step applied in MyLake C. The atmospheric equilibrium concentrations of CO 2 were calculated from the measured atmospheric CO 2 mixing ratios. The daily averages of the depth of the AML were estimated from the daily averaged temperature profiles as the depth at which water column temperature was within 0.25 • C of the temperature at 0.2 m as in Erkkilä et al. (2018). As in MyLake C, the approximation 20 U 10 /U 1.5 = 1.22 was used in the calculations. Following Mammarella et al. (2015), a value of 2 m −1 was used for the total attenuation coefficient of shortwave radiation K L in the calculation of Q eff .

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 25 new versions of the MyLake C Lake Kuivajärvi application, using the models for k by Cole and Caraco (1998)  The calibrations were performed against the daily averages of the automatic water column CO 2 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 2 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 2 fluxes. Even though the near-surface CO 2 concentration was the most significant factor considering air-water CO 2 exchange, deeper depths were included so that model behavior would remain reasonable also 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 2 , the 5 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 DOCrelated specific attenuation coefficient of photosynthetically active radiation β DOC , the maximal phytoplankton growth rate at 20 • C µ 20 , the phytoplankton death rate at 20 • C m 20 , the degradation rates of labile DOC k DOC,1 and semilabile DOC 10 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 15 included 1500 parameter sets. The medians of the final posterior distributions (Figs. S1-S4) were chosen as the calibrated parameters. They are presented, together with the default parameters, in Table 2. After the calibrations, additional goodnessof-fit metrics were calculated. The Nash-Sutcliffe efficiency (NS) 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 * > 0) or underestimation (B * < 0) of a state variable in the 20 simulation, whereas the normalized unbiased root-mean-square difference (RMSD * ) shows if the standard deviation of the simulated values is higher (RMSD * > 0) or smaller (RMSD * < 0) than that of the measurements (Los and Blaas, 2010).

Calculation of CO 2 budgets
After the calibrations, we calculated CO 2 budgets for the epilimnion of the lake during the periods of continuous summer stratification in 2013 and 2014 for each GEM. The epilimnion was defined as the layer in which water temperature was within 25 1 • C of surface temperature. The stratified period was defined to begin on the day of the formation of the thermocline after ice-off and to finish when the depth of the epilimnion (z epi ) reached the value of 7 m in the simulations. The exchange of CO 2 between the epilimnion and the atmosphere is balanced in MyLake C by (1) net external loading of CO 2 , (2) net epilimnetic CO 2 production, and (3) the release of CO 2 from deeper layers to the epilimnion. The net external loading equals the amount of terrestrially produced CO 2 entering the lake via stream inflow subtracted by the amount of CO 2 in lake outflow. The release 3 Results

Model calibration
Even though the differences between the formulations of the gas exchange models incorporated into MyLake C are rather notable, the resultant CO 2 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 5 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 2 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 2 concentration. The performance metrics for CO 2 concentration shown in the Supplement (Table S1)  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 2 concentration was successful also 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 25 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 , resisted wind-induced thermocline erosion during summer.

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 correspond- 30 ing 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 few occasions in early May and in October, whereas measurementbased 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 2 = 0.39-0.41, RMSE = 48.2-49.2 W m −2 , NS = 0.11-0.14, B * = −0.47 . . . − 0.46, n = 164). The differences in the simulated Q eff 's between GEMs, resulting mainly from different surface temperatures, were quite small. 5 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 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 10 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 2 = 0.35, RMSE = 0.11 m s −1 , NS = −3.2, B * = −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 15 1.589 × 10 −3 to 1.593 × 10 −3 between the GEMs.

CO 2 exchange
The differences between simulated gas transfer velocities for CO 2 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 include also the effect of thermal convection on gas exchange (Fig. 3, 20 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 's 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 25 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 2 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. 30 The simulated near-surface CO 2 concentrations were significantly too low during most of the study period in all GEMs except for CC, which yielded too high concentrations in autumn (Fig. 4, Table S3). The higher the simulated k and daily CO 2 efflux, the greater was the resulting decrease in near-surface CO 2 concentration. The decline of near-surface CO 2 concentration after ice-off was too rapid in all GEMs, especially in MI, the GEM with the highest k. Simulated near-surface CO 2 concentra-tion declined close to the atmospheric equilibrium concentration in all GEMs also in late summer because of insufficient gain of CO 2 in a shallow epilimnion developed under warm and calm conditions. The low simulated air-water CO 2 concentration gradients in May resulted in an underestimated air-water CO 2 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 . On the contrary, CC notably overestimated the corresponding calculated CO 2 flux in August 5 and September because of a high simulated near-surface CO 2 concentration. Also, the simulated CO 2 flux was slightly higher than the calculated flux in HE in August and September because of high epilimnetic net CO 2 production. The total simulated CO 2 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 2 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 10 and calculated fluxes in HE and MI.
The applied gas exchange models yielded notably different calculated monthly CO 2 effluxes ( Table 3). The CO 2 fluxes were calculated using the measured air-water CO 2 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 of those calculated with the other wind-based model CC. Days with a positive β, resulting in a lower k MI , occurred mainly in May and 15 October, and thus the difference between the CO 2 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 2 fluxes than the simplest model, CC. Nevertheless, the CO 2 fluxes calculated with MI were slightly higher than those calculated with HE. The CO 2 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 2 flux were dependent on the calculation interval. If the daily k had been calculated 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.

Lake CO 2 budgets
The simulated CO 2 budgets for the epilimnion of the lake during the periods of continuous summer stratification in 2013 and 2014 differed between GEMs as a response to different CO 2 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 rather well the depths estimated from the measured 5 daily temperature profiles (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 2 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 2 in the epilimnion decreased during the stratified period in 2013, whereas it increased slightly in 2014.
A higher net in-lake CO 2 production or a higher terrestrial CO 2 load was required to compensate for the higher CO 2 efflux in 10 GEMs that yielded higher values of k (Table 4) However, a high phytoplankton biomass did not imply high CO 2 consumption because of phosphorus limitation of phy-20 toplankton growth in the model and the resultant reduction of photosynthetic CO 2 consumption under high Chl a and low bioavailable phosphorus concentrations in the simulations. Instead, CO 2 fixation occurred at a steady rate and the total CO 2 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 2 fixation; however, also total net CO 2 production was 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 2 25 production despite low CO 2 fixation and a high k POC,sed in TE. Net CO 2 production was lowest in CC because of rather high total CO 2 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 2 input to the lake in GEMs with a high CO 2 efflux. The measured inflow CO 2 concentration was However, the total net external CO 2 loads to the lake over the stratification periods were slightly higher than the net external CO 2 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 5 in 2013 and 2014, respectively, the proportions being highest in CC and lowest in MI. The amount of CO 2 outflow was relatively large in CC because of the high epilimnetic CO 2 concentration; thus, the net external CO 2 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 2 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 10 in the CO 2 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 2 effluxes between GEMs.

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

25
The day-to-day performance of the simulation of epilimnetic CO 2 concentration was also partly determined by the simulated thermal stratification and epilimnetic volume. The simulations generally yielded too low a near-surface CO 2 concentration 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 2 concentration when the epilimnion became thicker, and vice versa, during the stratified period in 2013. Thermocline tilting-induced upwelling 30 and convection-induced entrainment transported more CO 2 -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 2 concentration. High solar radiation also enhances photosynthesis and thus increases the uptake of CO 2 (Provenzale et al., 2018). An overly deep simulated epilimnion resulted in enhanced CO 2 release from deeper layers and a higher total net CO 2 production in a larger epilimnetic volume, which were able to compensate for the CO 2 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 the gas exchange models are used either to calculate or to simulate daily estimates of k. The calculated 5 half-hour Q eff was generally directed into the lake on some occasions at daytime because of solar heating of the AML and always directed out of the lake at nighttime, and z AML often increased during nighttime and decreased under radiative heating of near-surface water at daytime. Boundary layer models and surface renewal models have been developed to describe shortterm dynamics of turbulence in a shallow AML, and thus they may not perform equally well in calculations with a daily time step.

10
The wind-based CC yielded the lowest and the surface renewal model TE the highest calculated air-water CO 2 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 2 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.

15
Because the simulations yielded significantly lower u * a compared to the values obtained through EC measurements (Fig. S8), the CO 2 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 20 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 25 and Prairie, 2013).

Comparison to EC CO 2 flux measurements
Estimates of air-water CO 2 fluxes obtained with the gas exchange models applied in our study have been compared with 30-min block-averaged EC CO 2 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's calculated with HE, CC, and MI with those obtained  Mammarella et al. (2015). In our study, the best agreement with simulated and calculated CO 2 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 2 fluxes, provided that the conclusions from the half-hour comparisons in 5 the above-mentioned studies can be extended to a daily scale.
The simulation results for the daily air-water CO 2 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 2 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 shorttime system malfunction may cause significant data loss during long study periods. Daily average or median EC CO 2 flux may 10 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 advection of CO 2 from the surrounding forest (Erkkilä et al., 2018). EC CO 2 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 2 concentration gradient due to the absence of photosynthesis as a CO 2 sink (Erkkilä et al., 2018). 15 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, however, may not represent the average conditions over the lake and hence may not yield correct estimates of whole-lake CO 2 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 . Thus, the values may not be representative for the whole lake. Wind speed and the resulting u * a 20 over lakes surrounded by forests are lower in sheltered near-shore 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 near-shore 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 25 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 30 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 2 concentration in small, shallow boreal lakes (Natchimuthu et al., 2017), which may result in further discrepancies in the estimates on whole-lake CO 2 flux obtained on the basis of gas exchange models or by using a vertical, horizontally integrated lake model.

Factors influencing the epilimnetic CO 2 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 2 concentration and air-water CO 2 flux, the different 5 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 2 dynamics.
Phytoplankton is an significant factor in the lake CO 2 budget and the main driver of the diurnal variation of CO 2 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 bacte-10 ria. 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 2 via degradation of phytoplankton-originated organic matter, as well as the release of bioavailable phosphorus in the epilimnion through 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 2 in the epilimnion was rather low in MI (Table 4)  A conspicuously high C DI,IN was needed to balance the high CO 2 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 2 , and the summertime inflow CO 2 concentrations may have been unnaturally high especially in TE. However, the use of 30 C DI,IN can be thought as the inclusion of the input of CO 2 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 2 concentration in groundwater in southern Finland is around 700-900 mmol m −3 (Lahermo et al., 1990), which is about tenfold higher than the estimated average inflow CO 2 concentration in Lake Kuivajärvi over the stratified period in 2013, 86 mmol m −3 , and well in line with the yearly average of groundwater CO 2 concentration near a boreal stream determined by Leith et al. (2015). Thus, also groundwater-derived CO 2 5 transport to the lake may affect the water column CO 2 concentration.
The effect of CO 2 inputs through minor inlets or groundwater may be supported by the fact that the simulated near-surface CO 2 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 groundwater level is generally relatively high (Fig. 4). The simulated epilimnetic CO 2 sinks were rather small at that time because net CO 2 consumption by phytoplankton was low in cool water 10 and because CO 2 efflux was relatively low because of a low air-water CO 2 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 2 concentration being approximately twice the simulated epilimnetic CO 2 concentration and the scaled inflow CO 2 concentrations and terrestrial CO 2 loads being even higher, the decline of epilimnetic CO 2 concentration was rapid in all GEMs. The high abundance of diatoms in Lake Kuivajärvi in spring may have resulted in 15 a supply of easily degradable organic matter, but net primary production also consumed CO 2 . Thus, substantial CO 2 loadings through surface runoff, minor inlet streams, or groundwater seepage could have been plausible additional sources of epilimnetic CO 2 in May, provided that the additional surface inflow was rich in CO 2 . 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 the discrepancy between the estimates of gain and loss of inorganic carbon was explained by a possible underestimation of the impact of groundwater inflow.

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 2 concentration and air-water CO 2 flux in the simulations. A higher CO 2 efflux would have required a higher gain of CO 2 in the lake through in-lake CO 2 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 25 the CO 2 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 2 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 30 made adaptable to higher CO 2 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 2 emissions and to study the development of lake ecosystems under changing conditions.

Conclusions
We studied the applicability of four gas exchange models with different complexity incorporated into a vertical physicobiogeochemical lake model MyLake C to the simulation of air-water CO 2 exchange and water column CO 2 concentration in a humic 5 boreal lake. The gas transfer velocities simulated using the simplest, wind-based gas exchange model by Cole and Caraco (1998), or CC, were best in accordance with the corresponding values calculated on the basis of direct in-lake measurements, whereas simulations with the other gas exchange models either overestimated (the models by Heiskanen et al. (2014) and MacIntyre et al. (2010)) or underestimated (the model by Tedford et al. (2014)) the respective calculated gas transfer velocities because of discrepancies in the simulation of wind stress or daily effective surface heat flux. 10 None of the applied gas exchange models resulted in a highly improved simulation performance regarding water column CO 2 concentration or air-water CO 2 flux. On the contrary, the more complex gas exchange models, which include both windinduced stress and heat-induced convection as the drivers of CO 2 exchange, yielded higher gas transfer velocities and thus higher CO 2 fluxes in the simulations, which resulted in difficulties in obtaining sufficient gain of CO 2 in the water column to balance the loss to the atmosphere. In addition, the model with a daily time step was not always able to simulate the changes 15 in near-surface CO 2 concentration and air-water CO 2 flux resulting from short-term physical processes, such as nighttime cooling or simultaneous surface heating and wind mixing. As a result, all the incorporated gas exchange models except for CC yielded notably too low summertime epilimnetic CO 2 concentrations in the simulations, which was also reflected by a significant underestimation of CO 2 fluxes compared to the corresponding fluxes calculated from the calculated gas transfer velocities and measured air-water CO 2 concentration gradients. The daily CO 2 fluxes simulated with CC were closest to the 20 corresponding calculated fluxes. The long and widely used CC was, however, shown to produce too low CO 2 flux estimates and its use was discouraged in an empirical gas exchange model intercomparison study by Erkkilä et al. (2018), whereas the more complex models yielded presumably more correct CO 2 fluxes. Nevertheless, it has to be noted that the comparison between the gas exchange models is more complex in our modeling study than in Erkkilä et al. (2018) because of the interplay between the simulated CO 2 flux and water column CO 2 concentration. 25 The present model application was not highly adaptable to increased CO 2 effluxes. The extent of in-lake production of CO 2 is largely related to model structure, process descriptions, and the estimation of parameter values, whereas the amount of external CO 2 inputs is governed by the quality of hydrological forcing data. Therefore, research on processes contributing to carbon cycling in boreal freshwaters and on the roles of different internal and external sources of CO 2 , such as groundwater, in lakes is sorely needed in order to enhance the predictive performance of model simulations. 30 The issues raised in our study concerning lacustrine carbon budgets can also be generalized to a larger scale. The application of advanced gas exchange models has been shown to lead to increased estimates of CO 2 emissions from boreal inland waters.
Thus, higher estimates of net terrestrial ecosystem production and carbon flux from land to inland waters are required to close the regional carbon budget. Also, the use of advanced, possibly more correct gas exchange models in the assessment of global gas efflux from freshwaters may result in higher estimates of the impact of freshwater ecosystems on global carbon cycling.
Code and data availability. The MATLAB model code for the MyLake C application presented in this study is freely available at https: //github.com/biogeochemistry/MyLake_C/tree/MyLake_C-gtsv. The automatic water column temperature and CO2 concentration data, the SMEAR meteorological data, and the manual measurement data presented in this study are available in the AVAA -Open research data 5 publishing platform (https://avaa.tdata.fi/web/smart/smear/download/). The metadata of these measurements are available via the ETSINservice (https://etsin.avointiede.fi/).
Author contributions. PK, TV, AO, and TH conceived the idea and designed the modeling study. PK developed the model application,