Lakes are disproportionately active sites in carbon cycling relative to their surface area (Cole et al., 2007), and with more accurate estimates of their numbers and surface area (Downing et al., 2006), emissions have been found to be significant relative to the terrestrial sources and sinks. The global estimate for the terrestrial uptake of CO2 is 3.0±0.9 Pg C yr−1 (Le Quéré et al., 2009), and the global evasion rate from lakes and reservoirs is now estimated at 0.32 Pg C yr−1 (Raymond et al., 2013). Many of the initial efforts demonstrating the contribution of inland waters to greenhouse gas evasion were based on measurements of fluxes with chambers on the water surface, by tracer studies, or by computing fluxes as the product of a gas transfer coefficient times the concentration gradient across the diffusive boundary layer at the air–water interface (Kling et al., 1992; Richey et al., 2002). In much of the recent effort to place emissions from lakes into the context of regional carbon budgets, computations were done using either a conservative value of the gas transfer coefficient or one based on wind speed (Cole and Caraco, 1998). Use of eddy covariance (EC) has led to upward predictions of gas evasion (Jonsson et al., 2008) relative to Cole and Caraco (1998) with some studies indicating that fluxes are enhanced under heating (McGillis et al., 2004) and others with cooling (MacIntyre et al., 2010). Similarly, more comprehensive models of the gas transfer coefficient indicate that improved accuracy may result from including other terms besides wind when estimating fluxes (Soloviev and Schlüssel, 1994; Soloviev et al., 2007).
The relative contribution of wind and buoyancy flux to turbulence at the air–water interface varies with lake size (Read et al., 2012). As small lakes are orders of magnitude more abundant than large lakes (Downing et al., 2006), further evaluation of their hydrodynamics as it applies to gas flux is required. For example, in northern latitudes many lakes are small, and in Finland >95% of the 190 000 lakes have an area <1 km2 and in total can cover locally up to 20% of the land area (Raatikainen and Kuusisto, 1990), and thus play an especially important role in regional carbon cycling. As turbulence mediates gas flux across the air–water interface (Zappa et al., 2007; MacIntyre et al., 2010; Vachon et al., 2010), the concern has been raised as to whether wind-speed-based parameterizations are applicable or need to be modified for small lakes (Read et al., 2012). By definition, at depths above the Monin-Obukhov length scale, L, which is the ratio of u*w3/0.4β, where u*w is the aquatic friction velocity, 0.4 is the von Karman constant and β the buoyancy flux (Imberger, 1985), wind-induced shear dominates turbulence production, whereas at depths below it, convection dominates. If that scaling is correct, wind shear would always be the dominant driver of emissions. Further, if turbulence is the dominant control on evasion at the air–water interface, a cubic dependency on wind speed is expected (Wanninkhof et al., 2009). Thus, with their numerical abundance in the landscape, understanding operative processes at the air–water interface in small lakes requires studies which quantify wind shear and buoyancy flux, and determines the relative contribution of each to the gas transfer coefficient.
Wind and cooling moderate gas transfer to the air–water interface. Their relative contribution has not been addressed for small lakes. The effect of wind on different strata can be predicted from the Wedderburn number, W, or its integral form, the Lake number, LN, which indicates the degree of thermocline tilting as a function of wind speed, the density gradient across the thermocline and basin morphometry (Imberger and Patterson, 1990). For high values of W, wind will only stir near surface waters. As W decreases below about 10, the thermocline begins to tilt. For W<3, mixing is anticipated near the margins; for W≤1, full upwelling is predicted plus considerable mixing throughout the thermocline and into the mixed layer (Monismith, 1985; MacIntyre et al., 1999; Horn et al., 2001; Boegman et al., 2005). The thermals produced during convection penetrate throughout the mixed layer and can cause entrainment of waters with higher concentrations of dissolved gases to the air–water interface (Crill et al., 1988; Eugster et al., 2003). Time series studies in small arctic lakes, with their reduced density difference across the thermocline, demonstrate that thermocline tilting and exchange between the upper mixed layer and more stratified layers can occur frequently, and mixed layer deepening and entrainment by convection are likely during cloudy conditions (MacIntyre et al., 2009). For boreal lakes, the density gradient across the thermocline will be larger in summer, but during fall cooling as it decreases, internal wave motions and convection could both contribute in bringing dissolved gases to the air–water interface and thus mediate increased fluxes.
EC systems, with their direct continuous measurements of flux (F), enable improved estimates of gas transfer coefficient, k (Eugster et al., 2003; Nordbo et al., 2011; Vesala et al., 2012). Using inverse procedures, the gas transfer coefficient k is computed as k=F/(caq−ceq), where caq is concentration of the dissolved gas in water and ceq is the concentration in the water in equilibrium with the atmosphere. The gas transfer coefficient is generally normalized to CO2 at 20°C and called k600. For the inverse procedure to be accurate in lakes, the footprint of the system must not reach onto the surrounding land (Eugster et al., 2003; Vesala et al., 2006), and the measurement of water concentration must be representative of the full footprint. In most cases, dissolved gases are only obtained at one location, although there are a number of studies which illustrate the spatial variability of dissolved gases (Roland et al., 2010; Rudorff et al., 2011; Abril et al., 2013). Spatial variability is expected if wind induces tilting of the thermocline, and some gases are mixed into the mixed layer either with the initial wind forcing or when the thermocline rebounds on cessation of the wind. Spatial variability would also result if the mixing from nocturnal cooling, which brings dissolved gases to the air–water interface (Crill et al., 1988), penetrated to different depths as would be expected if some parts of a lake basin were more exposed to wind and with concomitant higher evaporation rates. Lastly, processes such as differential cooling would bring gases produced nearshore to offshore sites where they could be mixed vertically by convection. Hence, evaluation of the accuracy of estimates of k600 from EC systems must consider whether spatial heterogeneity in gas concentrations has confounded the estimates.
Here we present results of a study using EC to quantify evasion of CO2 and the gas transfer coefficient in a small boreal lake during fall cooling. We developed an equation for the gas transfer coefficient based on boundary-layer theory, which took into account both wind speed and buoyancy flux. This equation is independent of the empiricism inherent in applications of the surface renewal model which currently rely on similarity scaling from ocean sites for estimates of turbulence (MacIntyre et al., 2001, 2010; Eugster et al., 2003; Read et al., 2012). We evaluated the new model and k from the surface renewal model and from regression equations in MacIntyre et al. (2010) relative to k from the EC system. With this approach, and computation of the Wedderburn and Lake numbers plus cooling, we were able to evaluate the accuracy of our modelling and assess whether gas concentrations were sometimes higher in the EC footprint due to quantifiable within lake motions. We evaluate the measured k600 during different conditions and compare it to the wind-based equation (Cole and Caraco, 1998) used in regional carbon budgets. With these combined approaches, we assess the physical processes which need to be included in modelling greenhouse gas evasion in small lakes.
The fieldwork was carried out in Finland in the humic Lake Kuivajärvi (61°50′ N, 24°17′ E) in Hyytiälä [Station for Measuring Forest Ecosystem–Atmosphere Relations (SMEAR) II]. The site is part of the ICOS (Integrated Carbon Observation System) measuring network, and the carbon and energy exchange between the forest and atmosphere at the site has been intensively studied (e.g. Kolari et al., 2009; Launiainen, 2010). The lake has an area of 63.8 ha and a catchment area dominated by a managed Scots pine (Pinus sylvestris L.) forest. Lake Kuivajärvi is a long and narrow north–south-oriented lake (Fig. 1), and with two distinct basins. The deeper south basin reaches a maximum depth of 13.2 m, making Kuivajärvi deeper than most lakes in Finland, since the mean value of the maximum depths of Finnish lakes is 8.2 m (Kortelainen et al., 2006). Lake Kuivajärvi is in the upper part of the Kokemäenjoki River basin, which drains to the Bothnian Sea in the northern part of the Baltic Sea.
The measurements were taken on a platform situated in the middle of the lake. The fetch is longest at the platform when the wind blows from 320° to 350°, approximately 1.8 km, or from 130° to 180°, approximately 0.8 km. From other directions, the fetch is only 0.4 km at maximum.
The measurements were carried out from 9 August to 30 November 2011; that is, they covered the end of the summer stratification period as well as the period of autumn overturn.
The EC measuring system on the platform was comprised of a Metek ultrasonic anemometer (USA-1; Metek GmbH, Elmshorn, Germany), which measured the three wind speed components and an enclosed path infrared gas analyzer LI-7200 (Li-Cor, Inc., Lincoln, NE, USA). The measuring height, zm, was 1.5 m above the water surface. The length of the polytetrafluoroethylene (PTFE) sampling tube was 1.8 m and the inner diameter was 8 mm. The flow rate was 12 L min−1, providing a turbulent flow. The tube was heated (3 W m−1) to avoid condensation. Sampling frequency was 10 Hz. The micrometeorological fluxes were calculated following Mammarella et al. (2009) and Aubinet et al. (2012). Several quality-screening criteria were applied to the EC data, including the flux stationarity test (Foken and Wichura, 1996), limits for skewness and kurtosis of scalar quantities and vertical wind velocity, as well as visual inspection of power spectra and cospectra. Much of the data during light winds and cooling did not pass the stationarity tests, zm/L>−0.55, and was not considered in the analysis. Therefore, most of the nightly data early in the period were missing. EC data post-processing and adopted quality criteria follow those used in Vesala et al. (2006) and Nordbo et al. (2011).
We obtained continuous measurements of water column concentrations of CO2 at depths of 0.2, 0.5, 1.0 and 3.0 m by using a similar measuring system as Hari et al. (2008), in which a continuous airstream was circulated with a diaphragm pump (KNF Neuberger Micro gas pump; KNF Neuberger AB, Stockholm, Sweden) in a closed loop. The system consisted of gas-impermeable tubing (stainless steel and Teflon), a CO2 analyzer (CARBOCAP® GMP343; Vaisala Oyj, Vantaa, Finland), semipermeable silicone rubber tubing with 1.5-mm wall thickness (Rotilabo 9572.1; Carl Roth GmbH & Co. KG, Karlsruhe, Germany) for gas collection and the pump. The semipermeable tubing was placed at the respective measuring depths. The gas collection tubing, CO2 analyzer and pump were connected with gas-impermeable tubing. The gas concentrations in the continuous airstream within the loop equilibrated with those in water around the semipermeable tubing; thus, the CO2 concentration in the water could be continuously measured in the gaseous phase. The silicone rubber was cleaned weekly during the open-water period by gently scrubbing the surface, and the rubber was entirely replaced every other month. The data were logged at 1-minute intervals. We additionally measured CO2 profiles at weekly intervals using the head space method (López Bellido et al., 2009).
For water temperature measurements, a chain of Pt100 temperature sensors was installed with sensors at depths of 0.2, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 8.0 and 10 m to continuously record water temperature. The sensors were new and calibrated at the factory. The accuracy of the sensors is 0.2°C and the response time 15 second, the data was recorded at 0.1°C resolution, the logging interval 5 second and the values at 0.2 m were used to represent surface water temperatures in the calculations. The incoming and outgoing longwave and shortwave radiation was measured with a CNR-1 net radiometer (Kipp & Zonen B.V., Delft, The Netherlands).
The gas transfer velocity was calculated from measurements as:
where is the EC CO2 flux, caq is the surface water CO2 concentration and ceq is the equilibrium CO2 concentration. Assuming that standard conditions prevail, the atmospheric pressure of CO2 can be calculated as follows:
where is the mole fraction of CO2 (given by the EC), and patm is the atmospheric pressure (=1 atm). The equilibrium concentration of CO2 (mol L−1) in the surface water is calculated using Henry's law:
where kH is the Henry's law constant at the surface water temperature. For ceq in g m−3, the results must be multiplied with the molar mass of CO2, thus resulting in the following final equation for equilibrium concentration:
The calculations for the bulk CO2 concentrations in water, caq, were made as with the equilibrium concentration in eqs. (2) and (4), and the mole fraction of CO2 was from the Vaisala probes. We used the data from the 0.2 m probe as the surface CO2 concentration, but we also computed expected CO2 concentration within the footprint taking into account the extent of thermocline deflection based on calculations of the Wedderburn and Lake numbers and our estimates of the depth of the actively mixing layer, and denote the gas transfer velocities obtained as ksite and kfp, where site stands for measurement site and fp for footprint. The Wedderburn number is a dimensionless ratio of the density stratification to the wind forcing and aspect ratio of the lake (Monismith, 1986):
where Δρ is the density difference between the epilimnion and the hypolimnion, ρ0 the mean density of the water column, g the acceleration due to gravity, h1 the depth of the epilimnion, L1 the length of the lake and u*w the friction velocity in water. A related index is the Lake number (Imberger and Patterson, 1990):
where S is the overall stability of the lake, zt the height of the centre of the metalimnion, H the depth of the lake, A the lake area and zg the height to the centre of the volume of the lake. LN, similarly to W, captures the influence of stratification, wind stress, and basin morphometry, and we use the indices identically in the following (MacIntyre et al., 1999). For W of order 1, CO2-rich metalimnetic waters would reach the surface. For partial upwelling with cooling either during or after the event, convection would mix water with higher CO2 concentrations into the upper mixed layer. For this purpose, we derived an equation which predicts the depth, zE, where E denotes for entrainment, at which the CO2 might reach the surface:
where h1/W is the maximum thermocline displacement (Monismith, 1986; Horn et al., 2001) and zAML the depth of the actively mixing layer. The zAML has been defined operationally as the depth in which the temperature is within 0.02°C of the surface temperature (MacIntyre et al., 2001; Rueda and MacIntyre, 2009). Due to the precision of our thermistors, we calculated zAML as the depth within the upper water column where the temperature was within 0.25°C of the surface temperature. In calculating corrected k, the CO2 concentration from zE was used as caq with a rough assumption that when the deeper CO2 mixed to the surface waters, the CO2 concentration would remain high for 7 hours. This time step was half the period of the dominant internal waves in the thermocline and approximately 2 hours longer than the transit time from the northern extent of upwelling to our measurement site assuming current speeds of 0.05 m s−1. Since Lake Kuivajärvi has two distinct basins, we used only the length of the southern basin in the W calculations, that is, L1=1.7 km, and the hypsographic curve for the southernmost basin in calculating LN .
The gas transfer velocity can be defined as:
where D is the diffusion coefficient for a given gas in water and z the mass boundary-layer thickness (Jähne and Haußecker, 1998). The thickness of the sublayer in which molecular diffusion exceeds turbulent diffusion can be defined as (Fairall et al., 2000):
where λ is an empirically determined coefficient, ν the kinematic viscosity of water, u* the friction velocity and Sc the Schmidt number.
The wind-induced water friction velocity and heat-induced water convective velocity were coupled as total water velocity by Jeffery et al. (2007):
where uref is the wind-induced water velocity at some reference depth, w* the waterside convective addition to water velocity and C2 an empirical coefficient. We wanted to provide a simple model, which could be implemented without the knowledge of, for example, surface roughness; therefore, in this study, we used a simple linear relationship between the wind speed and wind-induced water friction velocity:
where C1 is an empirical constant computed using a fixed drag coefficient and ratio of density of air and water, and U the measured wind speed.
When k is modelled as a function of wind speed, the water movement is assumed to be only wind-induced, as indicated by u* in eq. (9). When the effect of heat flux is also accounted for, the total water velocity, V, should be used instead of u*. Thus, the new equation for total k derived from eqs. (8)–(11) is (λ was set to unity):
The dimensionless constants C1 and C2 were fitted to the data with values 0.00015 and 0.07, respectively. The penetrative convective velocity w* was estimated according to Imberger (1985):
where β is the buoyancy flux, which was calculated as in Imberger (1985):
where α is the coefficient of the thermal expansion of water, Heff the effective heat flux and Cp the specific heat of water at constant pressure. Heff is the sum of latent and sensible heat flux, net longwave radiation and the portion of the shortwave radiation, which is trapped to the mixing layer (Imberger, 1985).
We applied the surface renewal model, k600=c1*(ɛv)0.25Sc−1/2, in which water friction velocity and buoyancy flux are used to estimate the rate of dissipation of turbulent kinetic energy, ɛ, and ν is kinematic viscosity, to estimate the gas transfer velocity (MacIntyre et al., 2010):
where the dimensionless constants c1, c2, c3 were fitted to the corrected k values and we obtained coefficients 0.50, 0.77 and 0.3, respectively, and the depth z was used as constant, 0.15 m. The u*w was calculated, accounting for atmospheric stability, as (Deacon, 1977):
where u*a is the atmospheric friction velocity, ρa the density of air and ρw the density of water.
We also used the model by Cole and Caraco (1998):
where U10 is the wind velocity at 10 m, derived from the measured wind velocity at 1.5 m, as U10=1.22U (Cole and Caraco, 1998).
The k value was also estimated from a regression model by MacIntyre et al. (2010), where the U10 has the units of m s−1, resulting in k in units of cm h−1:
All models were fit to the in situ temperature with the Schmidt number correction (Jähne et al., 1987).
As expected with fall cooling, air and surface water temperatures declined during the full study with air temperatures typically less than surface water temperatures (Fig. 2a). The concentration of CO2 in the surface water was initially below 1000 ppm (Fig. 2b), then fluctuated between 1000 and 1500 before reaching near steady values of 1500 ppm, and then declined again to 1000 ppm. These variations were related to lake physics, as will be discussed later. The concentration of CO2 in the atmosphere was near 360 ppm during daytime, rose to over 400 ppm during night in the first half of the measurements, and settled to ~390 ppm by the end of season. Occasionally, atmospheric CO2 concentrations of >500 ppm were measured at night. These occurred when winds rotated (Fig. 3f) and thus probably are an artefact of respiration from nearby forest. These values were discarded and k was not computed for such periods.
The components of the surface energy budget contributed in different ways to fall cooling. In general, the sensible heat flux was small, from −40 to 20 W m−2 (Fig. 2d). The sensible heat flux tended to be negative with some positive values when the lake was warming, for example, on days 238–240, 273 and 306. With positive values, we indicate a flux to the lake, while negative values signal flux to the atmosphere. The latent heat flux (Fig. 2d) varied between −300 and 0 W m−2 during the first month and became negligible when surface water temperatures were low. Maximum net shortwave radiation (Fig. 2e) was around 800 W m−2 during daytime at the beginning of the measurements, then steadily declined to <50 W m−2. The net longwave radiation (Fig. 2e) was always negative, occasionally up to −120 W m−2, but on average −32 W m−2. There were occasional rains (Fig. 2f), of which the most distinct event was the period of 2 weeks beginning on day 254. However, although it was raining every day during that time, the intensities were typically low, that is, <4 mm h−1.
The CO2 flux was negative all the time (Fig. 3a); that is, the lake was a source of CO2 to the atmosphere. In the first half of the measurements, the flux varied in general between −10 and 0 g CO2 m−2 day−1, with occasionally higher values. During the second half, the CO2 effluxes tapered off, in general <3 g CO2 m−2 day−1 after day 280. Overall, the CO2 efflux was <7 g CO2 m−2 day−1 about 90% of the time. The effective heat flux was largely negative during fall cooling with diel cycles of heating and cooling more common early in the study when solar radiation was higher (Fig. 2e, 3b). On average the effective heat flux was ca. −60 W m−2 in the first half, and −30 W m−2 during the second half of the study. The wind speed varied from 0.1 to 8.4 m s−1 (Fig. 3c). About 50% of the time the speed was <1.6 m s−1, 65% less than 2.0 m s−1, and <3.0 m s−1 about 85% of the time. There was significant channelling in the wind direction over the lake (Fig. 1); that is, 55% of the time the wind blew from the south (120°–180°) and 30% from the north (300°–360°), but there were periods when the wind rotated, for example, on days 222–226, 253–255 and 312–325 (Fig. 3d). Typical of these periods, wind speeds were also low. There were five periods of stronger winds when the wind speeds were close to or exceeded 6 m s−1. Three of the most distinctive periods were in the second half of the study, the first on days 286 and 287, second on 291–293 and the last on 308–310. The Wedderburn number was calculated for the period of seasonal stratification (Fig. 3e), with, as expected, low values when wind speeds were high. Before day 255, when wind speeds were >4 m s−1, W declined to <5, as is seen, for example, on days 221–222, 228, 233–234, 241 and 247, and on days 248–249 when U>5 m s−1, W was <2. Values of LN were slightly lower (data not shown) and supported the interpretation that internal wave induced mixing was likely whenever W<5. After day 255, as the lake continued to cool and the temperature gradient between epilimnion and hypolimnion reduced, W often dropped to <2 when winds rose >4 m s−1, that is, on days 255–256, 258–259, 262, 265–267, 273–277, 279 and 283.
The stability parameter zm/L (Fig. 3f), where L is the Monin-Obukhov length scale in the atmosphere, was typically between −2 (unstable conditions) and −0.01 (near-neutral conditions). There were some periods when the lake was heating and the zm/L was positive, being then in general from 0.2 (slightly stable conditions) to 0.001 (neutral conditions). High heat loss co-occurred with low winds speeds at night early in the study, when the atmosphere above the lake was unstably stratified, and zm/L was<−0.55. Most of these data at such times were discarded due to strict quality criteria. During periods of higher wind speed, near-neutral conditions (absolute value of zm/L<0.05) were observed.
In early August, there was a strong seasonal thermocline between 5 and 9 m depth with temperature difference across it of about 8°C. The temperature in the upper mixed layer was about 18°C (Fig. 4a). The seasonal stratification decreased as the lake cooled, and the mixed layer deepened. The CO2 concentrations in the mixing layer were around 800 ppm in the beginning of the measurements (Fig. 4b). There are four distinct periods of upwelling from the thermocline to the mixed layer on days 221–253. All of these events were related with W dropping to <5, and with concomitant cooling, the entrainment led to increased surface water CO2 concentrations, as seen on days 228–229, 233–235, 241 and 248–251.
Concentrations of CO2 in the surface water were variable. Increases were sometimes explained by upwelling, but at other times decreases occurred (e.g. day 258). One explanation could be dilution from rain fall. Alternatively, the decrease could be from advection. Wedderburn numbers dropped to 2 on day 256 and as the thermocline relaxes after such events, water from elsewhere in the basin can be advected into the measurement region contributing to variability in concentrations.
When the upper mixed layer deepened to over 8 m on day 266, the CO2 rich hypolimnetic waters mixed with epilimnetic waters, and CO2 concentrations increased to values >1400 ppm. After day 285 the whole lake water column was nearly isothermal, the water temperature was 7°C, and CO2 concentrations had decreased such that values were close to 600 ppm in the water column before ice on in early December.
We used two approaches to calculate the k using the EC data. The first was with the measured surface CO2 concentrations at the measurement site [eq. (1), Fig. 5a, dashed red line, ksite], and the second was based on the estimated surface CO2 concentrations in the EC footprint in which caq in eq. (1) are concentrations at depth zE estimated from the Wedderburn number (Fig. 5a, solid blue line, kfp), see section 2.3 for details. Both k were high in the first half of the measurements, ksite being on average 7.5 cm h−1 and kfp 6.4 cm h−1 and ranged from near zero to over 25 cm h−1. As a response to diurnal heating and nocturnal cooling, effective heat flux and atmospheric CO2 concentrations varied over diel cycles. During the stratified period, using CO2 concentration from depth zE reduced the number of times when k values were outside typical ranges (Banerjee and MacIntyre, 2004). We recorded during light winds 104 events with ksite when k>15, and only 49 similar events with kfp. When heat flux declined and variability decreased after day 293, the variation in k was reduced. Values of both ksite and kfp decreased towards the end of the measurements and were identical after day 285 when measured CO2 concentrations were near uniform. Two example periods of measured and modelled k are discussed in more detail later.
To determine conditions when buoyancy flux makes a considerable contribution to mixing relative to wind in the actively mixing layer, we contrasted convective and aquatic friction velocities (Fig. 5b). To determine whether buoyancy flux contributes appreciably to the gas transfer coefficient, we compared the wind shear term and the buoyancy term in our new model, kUw (Fig. 5c) and with the scaling in the surface renewal model, kSR (Fig. 5d). Early in the study, u*w and w* were similar in magnitude (Fig. 5b). Typically, at night during the first half of the measurements, w* was high, >0.004 m s−1, and u*w low, <0.002 m s−1, indicating that convective mixing would distribute dissolved gases throughout the mixed layer. In the daytime, the buoyancy flux was typically positive and the w* equalled zero by definition. When the heat flux decreased towards the end of the season, w* was <0.002 m s−1. During high wind events, the u*w reached or exceeded 0.01 m s−1, indicating wind was more effective in mixing as surface water temperatures cooled. The wind parameter in kUw was typically <0.5*10−6 m2 s−2, except when winds rose >5 m s−1, resulting in values >10−6 m2 s−2, as seen on days 249, 280, 286, 292 and 332, and all these events were also visible as peaks in the k data. The heat flux parameter in kUw followed the pattern of convective velocity (Fig. 5b), and when the parameter values were high, for example, >0.15*10−6 m2 s−2 during days 240–249, also the k tended to be high, around 10 cm h−1. Thus, with the scaling in kUw, the buoyancy flux term contributes to the gas transfer coefficient. The wind parameter in the surface renewal model responded to rising winds as seen for example, on days 249, 286 and 292 when winds of >7 m s−1 resulted in parameter values reaching or exceeding 2*10−11 m4 s−4. The heat flux parameter in kSR also followed the pattern observed in w*. However, it was always at least one order of magnitude smaller than the wind parameter, with only occasional exceptions very early in the study when both parameters were of equal magnitude. As a result of the different scaling, buoyancy flux had larger contribution in kUw than in kSR. The accuracy in the different weighting approaches between kUw and kSR are important for assessing the generality of the surface renewal model and will be discussed later.
Results from the models of k, that is kUw, kSR, kRA and kCC [eqs. (12, 15, 18 and 17)], are contrasted with ksite and kfp during stratified (Fig. 6) and unstratified (Fig. 7) conditions. Until day 244, the ksite and kfp were almost identical, and the values were close to the modelled k. Prior to 244, winds were southerly with a diel cycle with stronger winds in the day and weaker winds at night and the lake gained heat at daytime, sometimes even over 200 W m−2 and lost heat at night. On day 244, winds shifted to northerly and were steady for three days, U ~2 m s−1. Also, the lake was cooling continuously. Lake numbers dropped to a value of 2 on day 244, and high frequency non-linear internal waves occurred in the upper thermocline between days 244 and 246. These waves are indicative of wind speeds relative to stratification being high enough to cause thermocline tilting. Concomitantly, the actively mixing layer deepened (Fig. 4a), indicating entrainment of metalimnetic waters to the surface. When measured surface CO2 concentration was used in k calculations, many times very high k, k>15, were observed and the average value between days 244 and 247 was 10.7 cm h−1 with ksite (Fig. 6, dashed line). When we assume that entrainment had occurred at the EC footprint and CO2 concentrations had mixed and were equal to those at zE, almost all the very high k were gone, and kfp (Fig. 6, blue line) was similar to the values modelled, and the mean k dropped to 7.6 cm h−1, which was nearer the modelled average of 5.4 cm h−1 with kUw and kSR.
After this period, the winds started to rise, first reaching 5 m s−1 on day 248 and then one of the highest winds measured in this study, over 8 m s−1 on day 249. During both of these days, very high k were measured with ksite, but when the decrease in W and entrainment from convection was considered (Fig. 4b), kfp was more similar to modelled k than ksite. Overall, the models kUw, kSR and kRA performed very similarly during this example period, whereas kCC was about half of the other models’ estimates.
The second example (Fig. 7) is from late autumn (days 300–320), when the heat flux was already very low, around 30 W m−2 from day 300 to 313, after which it dropped to about −50 W m−2. During this period, with ksite and kfp identical, the measured k followed the patterns observed in wind speeds, for example, as seen on days 300–301 and 308–310 (Fig. 3c, 7 blue line) and when strong winds of >4 m s−1 prevailed, the measured k stayed high, >10 cm h−1.
During the second example, kUw performed quite similarly to kRA, and closely followed the measured k. The kSR quite constantly overestimated the measured k by 2 cm h−1 except in the event of high winds on days 309–310, when kSR was close to the measured k. The kCC was almost all the time <2 cm h−1 with the exceptions at strong winds, for example, on days 309 and 310 when it rose to 4 cm h−1. However, kCC was still underestimating the measured k by a factor of 2.
With respect to the study as a whole, the models kUw, kSR and kRA performed similarly and the comparison was better with kfp than with ksite. Virtually all of the very high k values with ksite were well above the modelled values, whereas when the CO2 concentration correction from entrainment was accounted for in kfp, many of the very high k decreased to the same range as the modelled k. This congruence supports our approach for estimating CO2 concentrations in the EC footprint. The regression model, kRA, and the surface renewal model had similar patterns to kUw and were more closely aligned when winds increased than during light winds with heat loss. Under such cases, kUw almost always provided the highest values, but regrettably these were also the occasions when most of the k data did not pass the quality tests and no reliable comparison could be made with the measured k. The kCC was quite steadily ~2 cm h−1, rising only when winds rose, and even during those conditions kCC was lower by a factor of two relative to other modelled and measured k. The average value for k was 4.2, 5.5, 4.2 and 2.2 cm h−1 with kUw, kSR, kRA and kCC, respectively, and 6.0 and 5.2 with ksite and kfp, respectively.
Comparisons of the modelled k against kfp (Fig. 8) showed similar performance for kUw (Fig. 8a), kSR (Fig. 8b) and kRA (Fig. 8c) and two-fold higher modelled estimates than with kCC (Fig. 8d). When kfp was <8 cm h−1, the relation to the modelled k (kUw, kSR, kRA) was close to linear, with greater linearity under heating conditions (data not shown). Once kfp exceeded 8 cm h−1, the modelled k values were underestimations. The kCC (Fig. 8d) underestimated kfp essentially all the time. The higher kfp relative to the models based on wind speed and buoyancy points either to the models not capturing an important aspect of the physics at the air–water interface or to continued underestimation of the CO2 concentrations within the EC footprint.
We obtained different relations for wind speed and k600fp, kfp after Schmidt number correction to CO2 at 20°C, during the stratified and unstratified period, when heat flux was negligible. During the stratified period (Fig. 9a), median values and the 25 and 75% values for k600fp showed a weak dependency on wind speed until winds reached 6 m s−1, with median k600fp from 5 to 9 cm h−1. This pattern was evident for both heating and cooling conditions, but the median values were ~5 cm h−1 under heating conditions and ~8 cm h−1 under cooling conditions (data not shown). We obtained the same results using Lake numbers with their slightly larger predicted thermocline tilt (data not shown). The trend is similar to McGillis et al. (2004), which predicted high continuous k even at low winds. In contrast, during the unstratified period (Fig. 9c), the data showed a power law dependence similar to McGillis et al. (2001), with median k600fp of 3 cm h−1 for U<3 m s−1. During seasonal stratification, k600fp had a linear relation to the effective heat flux during cooling conditions, kfp=−0.026Heff+6.0 (Fig. 9b), and no relation during heating. During autumn overturn, when the heat fluxes were quite small, no trend was visible between Heff and k600fp (Fig. 9d).
Evaluation of the relations between environmental parameters, the different models for k, and the kfp with the Pearson correlation (Table 1) showed good comparisons between models and a strong influence of wind in all seasons. For the full period, the correlations of the various models with kfp varied between 0.48 and 0.62 and were slightly higher than those with wind alone (Table 1). The correlations between the kfp and effective heat flux, w* and buoyancy flux were −0.47, 0.32 and −0.47, respectively, during cooling periods. All the models performed best during autumn overturn, when the correlation with kfp was ~0.6, and the correlation between kfp and wind speed was 0.61. The poorest correlations between the models and kfp were during seasonal stratification, when the environmental drivers, wind speed and heat flux, had low correlations with kfp, 0.29 and −0.32, respectively. Since there was a trend between heat flux and k600fp (Fig. 8b), the low correlations were related to high scatter in kfp, which was related to both the highly varying conditions within the lake ( Fig. 3e, 4) and EC random uncertainty, which is about 30% with 0.5-hour averaged CO2 fluxes. Again, the k data were mostly discarded during high heat loss and low winds, which confounds the correlation between kfp and heat loss expected from MacIntyre et al. (2010) and Read et al. (2012).
Our comparison of the relation between k600 and winds during the unstratified and stratified periods on Lake Kuivajärvi allows generalities which will apply to other lakes. During the unstratified period, the dependence on wind speeds occurred over a larger range of winds with the power law dependence greater than quadratic. These results point to turbulence controlling the gas transfer coefficient (Wanninkhof et al., 2009). The intercept for low winds is similar to the values in McGillis et al. (2001) and 50% higher than in Cole and Caraco (1998). During the stratified period, k600 did not depend on wind speed until winds reached approximately 6 m s−1. The median k600 values were of order 6–9 cm h−1 at wind <6 m s−1, similar to that observed by McGillis et al. (2004). The k600 values depend also on the magnitude of buoyancy flux under cooling and are independent of buoyancy flux under heating. The combination of rigorous quality control of our EC data, and our derivation of k which includes modification of gas concentrations induced by thermocline upwelling led to improved estimates of k600. Therefore, we conclude that the dependency of k600 on wind speed observed here and similar to observations in McGillis et al. (2001) is correct under unstratified conditions and stratified conditions with heating. In contrast, under cooling and for winds up to ~6 m s−1, median k600 values were of order 7 cm h−1. While these results are similar to observations of McGillis et al. (2004) in an oceanic study, it is uncertain whether the enhancement under cooling conditions was due to a stronger contribution of buoyancy flux than captured in our models or whether it was due to not fully capturing the spatial variability of CO2 concentrations in the EC footprint. Our measured mean gas transfer velocity for the whole study period was ksite=6.0 cm h−1, kfp=5.2 cm h−1, and kUw and kRA were ~4 cm h−1, that is, two times kCC of 2.2 cm h−1. These differences support recent work indicating results obtained using Cole and Caraco (1998) are underestimates (Rudorff et al., 2011), but see Schilder et al. (2013) and Vachon and Prairie (2013).
While wind was the dominant driver of the gas flux during overturn and at high winds, U>6 m s−1, winds were <6 m s−1 over 95% of the time. Our results indicate that there is uncertainty in estimates of k600 during cooling, and that continued effort to address this issue is required. Since low winds are a norm on small, sheltered lakes, and small lakes are orders of magnitude more abundant than large lakes (Downing et al., 2006), such effort is essential.
Our study indicates that Wedderburn numbers drop to values indicative of thermocline tilting and non-linear wave formation in small boreal lakes during fall cooling. Consequently, internal wave induced motions as well as convection from cooling deliver greenhouse gases as well as other dissolved substances to the upper mixed layer. Importantly, the resulting higher CO2 flux is captured by EC technique where the flux, and inherently the spatial variability of CO2 concentrations, from large areas is recorded. Thus, the EC technique, in contrast to chamber measurements, unless distributed over a lake surface, provides more accurate estimates of fluxes.
When lakes are stratified and internal waves moderate surface concentrations of greenhouse gases, accurate estimates of k depend on taking into account within lake physics when deriving k. When this effect was considered in kfp using the Wedderburn number and actively mixing layer depth, the measured mean k was reduced about 13% and many of the very high k, values >15 cm h−1, when neither wind nor heat loss indicated high k, decreased to values estimated with models which were independent from measured concentrations. Here we show, then, that the effects of wind and buoyancy flux at the surface mediate the gas transfer coefficient, and, when Wedderburn numbers drop to critical values, horizontal advection and vertical mixing from wind and buoyancy flux must be included in assessing fluxes of gases to the air–water interface. Consequently, focus should be put in developing more sophisticated, physically derived models and coupling these with biogeochemical models in order to develop generalities and for scaling up.
The mean CO2 flux during seasonal stratification was 1.15 µmol m−2 s−1, which was about four times the flux measured on another Finnish lake (Vesala et al., 2006; Huotari, 2011) and on a lake in Switzerland (Eugster et al., 2003), and almost five times the average flux from Alaskan lakes (Kling et al., 1992). During autumn overturn, the mean CO2 efflux was 0.56 µmol m−2 s−1, still higher than the efflux from the other lakes, indicating that boreal lakes surrounded by a managed forest might be an exceptionally high source of CO2 during the whole open-water period.
The results presented here provide solutions for improved estimates of k using EC approaches. Accounting for the spatial variability of surface CO2 concentrations allows for a more accurate derivation of k with the EC system. Here, we used a simple model of thermocline tilting due to decreases in the Wedderburn number below critical values. This approach significantly decreased the number and magnitude of high values of k. While other studies could similarly use this model, accuracy will be improved by measuring CO2 concentrations at multiple sites. We are uncertain as to the extent of error with spatially variable cooling, but it can also contribute to uncertainty and should be included in models. By consideration of the within-water column physics, estimates of k as a function of wind and buoyancy flux will be improved allowing for greater accuracy in estimates at global scales.
This study was funded through the Academy of Finland Center of Excellence program (project 1118615), Academy of Finland ICOS project (263149), EU ICOS project (211574), Research Foundation of the University of Helsinki, EU GHG-Europe project (244122), EU-project GHG-LAKE, DEFROST (Nordforsk) project, Academy of Finland (project 218094), and University of Helsinki research funds through the project Vesihiisi. Funding was provided by U.S. National Science Foundation Grants DEB 0919603 and ARC 1204267 to SM. We thank H. Miettinen for practical help during the field activities and the staff of the Hyytiälä Forestry Field Station for invaluable assistance and friendly working atmosphere.
Banerjee S , MacIntyre S , Grue J , Liu P. L. F , Pedersen G. K . The air–water interface: turbulence and scalar exchange . Advances in Coastal and Ocean Engineering . 2004 ; World Scientific, Singapore . 181 – 237 .
Cole J. J , Prairie Y. T , Caraco N. F , McDowell W. H , Tranvik L. J , co-authors . Plumbing the global carbon cycle: integrating inland waters into the terrestrial carbon budget . Ecosystems . 2007 ; 171 – 184 .
Downing J. A , Prairie Y. T , Cole J. J , Duarte C. M , Tranvik L. J , co-authors . The global abundance and size distribution of lakes, ponds, and impoundments . Limonol. Oceanogr . 2006 ; 51 : 2388 – 2397 .
Eugster W , Kling G , Jonas T , McFadden J. P , Wüest A , co-authors . CO2 exchange between air and water in an Arctic Alaskan and midlatitude Swiss lake: importance of convective mixing . J. Geophys. Res . 2003 ; 108 : 1 – 12 .
Hari P , Pumpanen J , Huotari J , Kolari P , Grace J , co-authors . High-frequency measurements of productivity of planktonic algae using rugged nondispersive infrared carbon dioxide probes . Limnol. Oceanogr. Methods . 2008 ; 6 : 347 – 354 .
Kortelainen P , Rantakari M , Huttunen J. T , Mattsson T , Alm J , co-authors . Sediment respiration and lake trophic state are important predictors of large CO2 evasion from small boreal lakes . Glob. Chang. Biol . 2006 ; 12 : 1554 – 1567 .
MacIntyre S , Eugster W , Kling G. W , Donelan M. A , Drennan W. M , Saltzman E. S , Wanninkhof R . The critical importance of buoyancy flux for gas flux across the air–water interface . Gas Transfer at Water Surfaces, Geophysical Monograph Series . 2001 ; Washington, DC : American Geophysical Union . 135 – 139 .
Mammarella I , Launiainen S , Gronholm T , Keronen P , Pumpanen J , co-authors . Relative humidity effect on the high-frequency attenuation of water vapor flux measured by a closed-path Eddy covariance system . J. Atmos. Ocean. Technol . 2009 ; 26 : 1856 – 1866 .
Nordbo A , Launiainen S , Mammarella I , Leppäranta M , Huotari J , co-authors . Long-term energy flux measurements and energy balance over a small boreal lake using eddy covariance technique . J. Geophys. Res . 2011 ; 116 : 1 – 17 .
Read J. S , Hamilton D. P , Desai A. R , Rose K. C , MacIntyre S , co-authors . Lake-size dependency of wind shear and convection as controls on gas exchange . Geophys. Res. Lett . 2012 ; 39 : L09405 .
Richey J. E , Melack J. M , Aufdenkampe A. K , Ballester V. M , Hess L. L . Outgassing from Amazonian rivers and wetlands as a large tropical source of atmospheric CO2 . Nature . 2002 ; 416 : 617 – 620 .
Roland F , Vidal L. O , Pacheco F. S , Barros N. O , Assireu A , co-authors . Variability of carbon dioxide flux from tropical (Cerrado) hydroelectric reservoirs . Aquat. Sci . 2010 ; 72 ( 3 ): 283 – 293 .
Rudorff C. M , Melack J. M , MacIntyre S , Barbosa C. C. F , Novo E. M. L. M . Seasonal and spatial variability of CO2 emission from a large floodplain lake in the lower Amazon . J. Geophys. Res . 2011 ; 116 : G04007 .
Schilder J , Bastviken D , van Hardenbroek M , Kankaala P , Rinta P , co-authors . Spatial heterogeneity and lake morphology affect diffusive greenhouse gas emission estimates of lakes . Geophys. Res. Lett . 2013 ; 40 : 1 – 5 .
Soloviev A. V , Donelan M , Graber H , Haus B , Schlüssel P . An approach to estimation of near-surface turbulence and CO2 transfer velocity from remote sensing data . J. Marine Syst . 2007 ; 66 : 182 – 194 .
Vachon D , Prairie Y. T , Cole J. J . The relationship between near-surface turbulence and gas transfer velocity in freshwater systems and its implications for floating chamber measurements of gas exchange . Limnol. Oceanogr . 2010 ; 55 ( 4 ): 1723 – 1732 .
Vesala T , Eugster W , Ojala A , Aubinet M , Vesala T , Papale D . Eddy covariance measurements over lakes . Eddy Covariance. A Practical Guide to Measurement and Data Analysis . 2012 ; Dordrecht, The Netherlands : Springer . 365 – 376 .
Vesala T , Huotari J , Rannik Ü , Suni T , Smolander S , co-authors . Eddy covariance measurements of carbon exchange and latent and sensible heat fluxes over a boreal lake for a full open-water period . J. Geophys. Res . 2006 ; D11101 .
Zappa J. C , McGillis W. R , Raymond P. A , Edson J. B , Hintsa J , co-authors . Environmental turbulent mixing controls on air–water gas exchange in marine and aquatic systems . Geophys. Res. Lett . 2007 ; 34 : L10601 .