Lake‐Atmosphere Heat Flux Dynamics of a Thermokarst Lake in Arctic Siberia

We conducted eddy covariance measurements from April to August 2014 on a Siberian thermokarst lake. The study site is located in the Lena River Delta and characterized as a floating ice lake. Heat fluxes differed in magnitudes, directions and temporal patterns depending on the lake surface conditions (“frozen” ice cover, ice cover melt, and open water). Significant heat release during frozen ice cover conditions highlighted the importance of lakes for the landscape heat budget and water balance. The energy balance was nearly closed during the open water period and highlighted the impact of melting energy on its closure during the ice cover period. Sensible and latent heat dynamics were driven by temperature and water vapor gradients scaled by wind speed, respectively. We calculated bulk aerodynamics transfer coefficients and evaluated the performance of the derived in situ and three independent heat flux parameterization schemes. We found that bulk transfer models perform moderately to poorly for the different lake surface conditions. During the open water period small‐scale temporal variability could not be represented by the models, particularly in case of latent heat flux. The model results were less sensitive to the specific model type than to the accuracy of the surface water temperature measurement, which is dependent on a well‐thought‐out measurement design. Our study stresses considerations that are crucial for similar campaigns in the future, in order to face the measurement challenges encountered on arctic lakes especially during the ice cover period.


Introduction
The Arctic is a key player for global atmospheric circulations and comprises potential tipping elements in the context of climate change (Duarte et al., 2012;Lenton et al., 2008;Wadhams, 2012). Due to polar amplification it experiences global warming at a much faster rate than the global mean (Arctic Monitoring and Assessment Programme, 2015; Duarte et al., 2012;Intergovernmental Panel on Climate Change, 2013;Jeffries & Richter-Menge, 2015), initiating feedback mechanisms such as increased carbon release due to permafrost thaw and degradation (Schuur et al., 2008(Schuur et al., , 2015Zimov et al., 2006). Thermokarst and thermal erosion further allow profound changes in the landscape including subsidence and drainage alteration. Consequently, the size and number of "migratory" thermokarst lakes are changing (Zimov et al., 2006), with great spatial variations depending on the ground-ice content Jorgenson et al., 2006;Morgenstern et al., 2011;Smith et al., 2005). Thermokarst lakes are of high importance to the atmosphere-subsurface thermal exchange in tundra environments. Their heat storage capacity can alter the heat budget of surrounding permafrost areas by transferring heat into the permafrost, initiating further permafrost thaw (Arp et al., 2016;Boike et al., 2015;Langer et al., 2016). Taliks develop underneath floating ice lakes with lake bed temperature above freezing and extend for a multiple of a lake's depth into the ground. The above-average warming of the Arctic and related reduced ice growth has already resulted in a shift in shallow thermokarst lakes from bedfast to floating regime in some regions, creating patchily thawed surfaces in continuous permafrost lowlands (Arp et al., 2012(Arp et al., , 2016. Lakes generally differ from vegetated surfaces in terms of their energy flux dynamics, with higher absorption of solar radiation during the open water period and a large heat storage capacity and associated inertia (Mammarella et al., 2015;Nordbo et al., 2011). Lake heat storage is regulated at the air-water interface by turbulent heat fluxes and long-wave radiation (Nordbo et al., 2011). Turbulent heat fluxes further influence the water balance and alter the strength and duration of thermal stratification and mixing, highlighting them as important features in the functioning of a lake (Verburg & Antenucci, 2010;Woolway et al., 2015).
A thorough knowledge and adequate parameterization of land-atmosphere interactions including the exchange of heat and moisture is highly relevant for reliable numerical weather prediction and climate modeling (Beyrich et al., 2006). In this context, lakes are gaining progressively more attention because of the increasing spatial resolution of the implemented land surface models. Depending on their proportional coverage within the landscape, they affect climate and weather through changes in surface energy budgets at distinct spatiotemporal scales (Dutra et al., 2010;Jeffries & Richter-Menge, 2015;Martynov et al., 2012). However, due to insufficient data, lake-atmosphere interactions are not yet satisfyingly understood and implemented in models (Martynov et al., 2012;Mironov et al., 2010). Remoteness and associated challenges in terms of infrastructure, maintenance, and safety result in particularly high model uncertainties for the Arctic, where lakes and ponds capture a substantial portion of the land surface area (Muster et al., 2013(Muster et al., , 2017Rouse et al., 2005). The ice cover period is particularly understudied, though it plays a key role by dominating the lake-atmosphere exchange dynamics for 60-80% of the year (Arp et al., 2011;Brown & Duguay, 2010;Prowse et al., 2015). Furthermore, existing studies on surface heat fluxes of arctic lakes are biased toward the North American continent (e.g., Arp et al., 2015;Potter, 2011). A few studies do however exist in arctic Scandinavia (Leppäranta et al., 2016) and Russia (Boike et al., 2015).
The eddy covariance (EC) technique has been employed for direct measurements of turbulent heat fluxes over lakes for about 30 years (Blanken et al., 2000;Stannard & Rosenberry, 1991;LITFASS-2003experiment, see Beyrich et al., 2006; see Aubinet et al., 2012 andMammarella et al., 2015 for an overview of recent studies; Erkkilä et al., 2018). However, there is a particular shortage of sites in the Arctic. Instead, less reliable indirect flux estimation methods are applied, which are based on more frequently measured variables such as water temperature (T w ) and wind speed (U). However, in the Arctic even these standard measurements are very scarce, especially in northern Siberia. Satellite observations can be utilized for the derivation of the surface heat fluxes but only for large lakes (Moukomla & Blanken, 2017).
One of the most common models used for indirect flux estimations is the bulk aerodynamic transfer model (flux-gradient method). Only a few studies employing this method exist for arctic thermokarst lakes (Huang et al., 2016;Potter, 2011) and commonly used transfer coefficients are not yet sufficiently validated for this lake type as respective measurement data sets are missing.
Bulk aerodynamic transfer models were initially developed in ocean research and applied either without modifications or only tentatively modified for lakes and reservoirs (Woolway et al., 2015). However, physical conditions over oceans and lakes can differ remarkably depending on the size of the lake (Verburg & Antenucci, 2010;Woolway et al., 2015). All parameterization schemes use the same basic bulk equations, where the turbulent heat fluxes are assumed to be linearly proportional to transfer coefficients but vary in the derivation of these coefficients (Brunke et al., 2003). Complex approaches derive the transfer coefficients with different roughness length parameterizations and treatments of atmospheric stability, free convective conditions, and surface layer gustiness (Verburg & Antenucci, 2010;Woolway et al., 2015).
In this study we aim to contribute to the understanding of the processes and magnitudes of the sensible (H) and latent (LE) heat exchange between an arctic thermokarst lake and the atmosphere by means of EC measurements. We will therefore investigate how the flux dynamics differ between contrasting lake surface conditions. We will further study how distinct components contribute to the energy balance (EB) of the lake and whether we can close this balance. We apply bulk transfer relationships for the indirect parameterization of the heat fluxes and evaluate the performance of available schemes for the different lake surface conditions. Our results aim to contribute to a better representation of the physical processes and dynamics of arctic lakes within regional and global modeling in order to reduce their uncertainties.

Study Site
Our study site, an arctic Siberian thermokarst (thaw) lake, is located on Kurungnakh Island, in the southerncentral part of the Lena River Delta at 72°17.9 0 N, 126°10.4 0 E (altitude 44 m; Figure 1). The fan-shaped river delta is the largest delta in the Russian Arctic, covering 25,000 km 2 . It comprises more than 1,500 islands and around 60,000 lakes on three distinct river terraces of different age and multiple floodplain levels. The Lena River Delta is located in the zone of continuous permafrost. On Kurungnakh, large thermokarst lakes and basins are essential components of the landscape (Morgenstern et al., 2011). The thermokarst origin of lakes is common in the lowland tundra permafrost areas of Northeast Siberia, characterized by fine-grained, organic-rich Yedoma-type ice complex deposits (Schirrmeister et al., 2013), similar to the arctic lowlands in central and eastern Siberia, the interior and north of Alaska, and northwestern Canada .
The study site is characterized by a dry continental arctic climate with very low temperatures and precipitation. Boike et al. (2013) reported an annual mean air temperature of À12.5°C and July as the warmest month (with a mean temperature of 10.1°C) for 1998-2011 on nearby Samoylov Island. About 45% of the precipitation is rain during the growing season (Boike et al., 2008). The 24-and 0-hr daylight periods last from 7 May to 7 August and from 15 November to 28 January, respectively.
The lake studied is one of 549 lakes, which cover about 7.5% of Kurungnakh Island (Morgenstern et al., 2011). It has an area of 1.21 km 2 and volume of approximately 3.8 Mio·m 3 with mean and maximum lake depths of 3.1 and 6.5 m, respectively. The lake is ice covered for 8-9 months each year. It is characterized as a floatingice and polymictic lake, with only a few days of summer stratification. Ice-on and ice-off were observed around end of September/beginning of October and June/July, respectively (Boike et al., 2015). A complete snow cover of the landscape is often only for a short duration, due to generally low snowfall in this region and strong winds resulting in local snow drifts and hardened wind crusts.

Field Measurements
In April 2014 we positioned a floatable measurement platform on the ice cover in the center of the lake. The quadratic 9-m 2 platform consisted of nine black plastic pontoons (Table S1 in the supporting information). Power was supplied by four solar panels on the corners of the raft. We equipped the platform with EC instrumentation (CSAT3, Campbell Scientific Inc.; LI-7200, LI-COR Biogeosciences) mounted on a tower and measured turbulent fluxes of momentum, heat, and H 2 O. We recorded raw turbulence and dry mole fraction data at 20 Hz in half-hourly files. The average measurement height during the study period was 2.39 m.
In addition to EC instrumentation, sensors for basic atmospheric variables (air temperature, T a , and relative humidity, RH) and platform position, as well as a time-lapse camera, were installed on the tower. Snow depth and radiation (incoming and outgoing long-and short-wave radiation LW in , LW out , SW in , and SW out ) sensors were installed on an arm branching out from the tower. A water temperature (T w ) profile was measured with a thermistor chain, which was attached to one of the platform corners, comprising eight T w sensors (intended measurement depths: water surface, 0.5, 1, 1.5, 2, 3, 4, and 5 m). Another T w profile was measured from an independent buoy (intended measurement depths: water surface, 0.5, 1, 1.5, 2, and 2.5 m, lake bottom). The buoy was initially placed in the northern central part of the lake and found on 9 August 2014 on the lake shore. Both profiles were measured every half hour.
Our measurement period lasted from 23 April to 16 August 2014, at which point the platform was dismounted. We defined three subperiods: the "frozen" ice cover period (23 April to 21 May and 27-29 May 2014), the "melting" ice cover period (including ice breakup; 21-26 May and 30 May to 24 June 2014), and the open water period (24 June to 16 August 2014) starting with ice-off. A 3-day cold period at the end of May (27-29 May 2014) was classified as the first subperiod as the ice surface was frozen instead of melting (section 3.1). During ice cover melt, the platform insulated the ice cover below, resulting in the elevation of the platform relative to the adjacent melting ice, until it slid into the water on 19 June 2014. Since then, the platform floated and rotated depending on the wind direction (WD; Figure S2b). However, the platform was already drifting and rotating together with the whole ice cover as soon as the nearshore shallow areas were ice free. The data set is shown in local time, which was 10 hr ahead of UTC time during the study period.

Heat Flux Computation
We calculated half-hourly fluxes of H and LE based on the block-averaged covariances between the respective scalar s and the vertical wind velocity w using the software package EddyPro® (Version 6.1) [Computer software] (2015) (LI-COR, Nebraska, USA; Text S1 for details on the EC processing; (Charuchittipan et al., 2014;Eugster et al., 2003;Foken, 2008;Fratini et al., 2012;Horst & Lenschow, 2009;Kaimal & Gaynor, 1991;Mammarella et al., 2009Mammarella et al., , 2016Moncrieff et al., 2004;Nordbo et al., 2014;Peltola et al., 2014;Ross & Grant, 2015;Siebicke et al., 2012;van Dijk et al., 2004;Vickers & Mahrt, 1997;Wilczak et al., 2001)). Quality control and filtering resulted in a final data coverage of 91% and 82% for H and LE, respectively. We calculated evaporation E for each half hour of the open water period in accordance with Verburg and Antenucci (2010). The density of water and latent heat of vaporization of water (L v ) were both calculated as functions of the near-surface (bulk) water temperature measured at the raft ( T w surf ), in accordance with Andreas (2005) and Verburg and Antenucci (2010), respectively.

Footprint Estimation
For footprint analysis we employed the two-dimensional parameterization for Flux Footprint Prediction (FFP), which is described by Kljun et al. (2015). We used the openly available R code v. 1.3 (http://footprint.kljun.net; Text S2 for the estimation of the planetary boundary layer height). As the platform changed its position during the study period, we calculated the cumulative footprint climatology for each of the main locations ( Figure S2c). The lake contributed at least 90% to the cumulative flux, except for the 3-week period, when the platform was located close to the eastern shore (10-30 July 2014), where it contributed at least 80%.

Water Temperature Profile and Heat Storage Change of the Lake
Based on the two thermistor chains, we documented the thermal characteristics of the lake throughout the study period. As the raft was progressively elevated from the surface during ice cover melt, the positions of the water temperature sensors did not exactly represent the intended depths until the platform was floating (19 June 2014), exposing, for example, the upper sensor to the air during frozen and melting ice cover. For the heat storage change calculation, we used the water temperature measurement (raft) at 0.5-m depth (T w0:5m ) as T w surf until the platform was floating and replaced these data with the water temperature measurements of the buoy at 0.5-m depth for a 2-day period during ice cover melt, when the raft reached its most highly elevated position, exposing this sensor.
Before calculating the lake heat storage change, we first created 3-hr-long running means of the raft T w measurements to lower their noise level and avoid the associated random errors (Nordbo et al., 2011). We then calculated the depth-weighted time derivative of the ice and water column temperature and the heat storage change of the lake. Although the T w measurements partly represented the ice temperature instead of the water temperature, we decided to calculate ΔQ exclusively based on the density and specific heat of water as T w was above 0°C. The density of water and specific heat of water were calculated in accordance with Andreas (2005) as a function of the depth-weighted average temperature T w (Gill, 1982).

Energy Balance Closure
We investigated the EB at the lake surface in order to evaluate the reliability of the EC heat flux measurements following Nordbo et al. (2011) and including net radiation (R n ), lake heat storage change (ice and water column; ΔQ), sensible heat flux (H) and latent heat flux (LE; all in W/m 2 ). We expected other components to be negligible, at least during the open water period, and respective field measurements were not conducted. The turbulent fluxes are positive when directed upward from the lake surface to the atmosphere, whereas R n and ΔQ are defined positive when directed toward the lake surface.
We found the net outgoing radiation to be underestimated due to an insufficient distance between the radiation sensor and the platform. We expected the platform to have reduced SW out and increased LW out due to its black color. In order to correct the SW out measurements, we estimated the source area contribution of the platform within the field of view of the pyranometer and subtracted this contribution from the measured albedo assuming an albedo of 0.02 for the platform (Text S3). Compared to the bias in SW out , the one in LW out is impossible to properly estimate based on the available data from our campaign.
In order to quantify the EB we calculated the ratio between the turbulent fluxes and the available energy (EB closure [EBC], in %) and the absolute EB residual Res (W/m 2 ) based on the respective cumulated EB components of an average day of each subperiod (mean daily course): Res ¼ R n À ΔQ À H À LE: We used all available half-hourly values, even if all the required EB components were not available. The EB is closed with EBC = 100% and Res = 0 W/m 2 . Compared to an EB calculation on half-hourly basis, this method has the advantage to compensate the intradaily dynamics of ΔQ, avoiding phase lags in the energy storage (Mammarella et al., 2015).
We assumed that a small component of available energy is used to melt the ice during the frozen ice cover period and a large component during the melting ice cover period, resulting in a small and large Res, respectively. To examine this, we roughly back-calculated the ice thickness at the beginning of the measurement period, assuming that Res corresponds entirely to the melting energy (Text S4). This latent heat of fusion represents the energy taken for phase transition from ice to water, provided especially by solar radiation but also heat flux from the water column to the ice and the turbulent heat fluxes.

Controls of Heat Flux Dynamics
We tested the correlations of H and LE with the gradient of the surface water temperature (T 0 ) and T a and the vapor pressure gradient multiplied by U, respectively (UΔT and UΔe). We calculated Δe according to Verburg and Antenucci (2010) and accounted for respective ice cover and open water conditions during the periods by separately deriving the saturated vapor pressure for ice and water, both according to World Meteorological Organization (2008). Only EC flux measurements of best quality (QC = 0; Mauder & Foken, 2004) were used for this analysis. The resulting linear regressions supplied essential input for the in situ bulk aerodynamic transfer model in section 2.8.1. We investigated all three subperiods separately as we expected them to be characterized by distinct heat flux dynamics due to the differing surface conditions. For comparison, we further derived the regressions when considering the whole study period.
Usually, skin temperature T 0s (calculated from measured LW out based on the Stefan-Boltzmann law, infrared emissivity of 0.98) has to be used as T 0 for this analysis as the heat exchange with the atmosphere takes place at a thin surface layer. For open water conditions, this layer is described as "cool skin" (Fairall et al., 1996). The use of T 0s was particularly important for the frozen ice cover period due to the significant linear vertical temperature gradient in the ice cover. However, the platform-induced overestimation of LW out (section 2.6) was directly reflected in an overestimated T 0s , which led to a bias in ΔT and Δe. As a consequence, the heat flux dynamics during the melting ice cover and open water period could not be reflected sufficiently. During these two subperiods we instead used the water temperature measured closest to the surface T 0m (T w0:5m during the melting ice cover and T w surf during the open water period) as the ice and water column was essentially isothermal (exceptions, e.g., at nighttime during ice cover melt). The temperature of melting ice and intermingled water is about 0°C.

Heat Flux Parameterization With Bulk Aerodynamic Transfer Models
We estimated the heat fluxes based on bulk aerodynamic transfer models on a half-hourly basis, comparing one in situ and three independent parameterization schemes, whereby the first is developed based on the measured fluxes and the latter are independent of measured fluxes. We used the same combination of T 0 and T 0s as in the calculation of UΔT and UΔe. The standard bulk transfer formulas are defined as where ρ a is the air density (calculated according to Verburg & Antenucci, 2010), C a is the specific heat capacity of air (calculated according to Andreas, 2005), C H and C E are transfer coefficients and L s, v is the latent heat of sublimation (for ice-covered conditions) or vaporization (for open water conditions). L s was calculated according to Andreas (2005) as function of T 0 .

In Situ Model
The bulk transfer coefficients of heat and water vapor C H, E are described by the slopes b H, E of the respective fitting curves for H and UΔT, and LE and UΔe over all stabilities and U and are estimated as The estimated transfer coefficients were applied in equations (3) and (4). We derived heat fluxes for the three subperiods separately and combined them into one time series . For comparison, we also parameterized the heat fluxes for the full study period (IN-SITU v2). As the coefficients were calculated for our specific measurement height and are thus site specific, we additionally estimated the coefficients in the reference height of 10 m (C H10 and C E10 ) following Xiao et al. (2013) from linear regression involving U at 10-m height (U 10 ). For an estimation of the roughness length z 0 , which was needed for the standard logarithmic wind law (neutral stability assumption) in order to estimate U 10 , we applied the flux gradient relations for momentum with stability classes as described in Woolway et al. (2015). The drag coefficient for momentum (C D ) was estimated as slope of the regression of the square of U versus the square of friction velocity (u * ).

Independent Models
The independent parameterization schemes applied in this study comprise two approaches for open water conditions, that is, lakes/reservoirs (LAKE1) and sea (SEA), and one specific approach for sea-ice (SEA-ICE1).
In comparison to the in situ parameterization approach, they consider the distinct atmospheric stabilities, with the transfer coefficients as functions of measurement height and atmospheric stability. The models are typically fed by T a , T 0 , RH, air pressure, U, SW in , and measurement height. Where necessary, we modified the model codes to accept flexible measurement heights. 2.8.2.1. Lake Heat Flux Analyzer (LAKE1) The Lake Heat Flux Analyzer is the GLEON (Global Lake Ecological Observatory Network) scheme for the calculation of surface energy fluxes in lakes (Woolway et al., 2015). It was developed for the analysis of highfrequency data from instrumented lake buoys deployed within Global Lake Ecological Observatory Network to facilitate robust comparisons between lakes due to standardized and consistent processing.  (Andreas et al., 2010). We applied this routine only for the frozen ice cover period, with the "winter" setting and an ice concentration (ratio of ice and water areas) of 100%, to test if this parameterization scheme developed particularly for ice-covered conditions provides better results than other schemes (model code sources are provided in Text S5).
We further tested three additional independent models for comparison (Text S5 for details).

Meteorological Conditions and Heat Flux Dynamics
At the beginning of our study period the lake was ice and partly snow covered. The ice cover was about 1.7 m thick (measured close to the platform). However, the T w profile derived from the thermistor chain attached to the raft indicated an ice thickness of about 2 m (Figure 2f). We observed intense snow drift during the frozen ice cover period due to high U (mean 4.9 m/s, max. 16.9 m/s; Figures 2c and S2a) and three events of fresh snow in April and May, as well as two events of rime in June. The dynamics of R n and surface albedo (Figure 2e) reflected the specific conditions of the lake surface during the ice cover. RH ranged between 32% and 100% with a mean of 85%; u * was on average 0.26 m/s without significant seasonal patterns but a very slight increase in the afternoon.
The heat fluxes differed between the subperiods and reflected changing stable and unstable conditions (Figures 2a and 2c). The frozen ice cover period was characterized by very low T a down to À16.6°C at the end of April (Figure 2d), remaining most often below T 0 and resulting in mean H of 10.8 W/m 2 ( Figure 2a; Table 1). With a mean flux of 18.6 W/m 2 LE was almost entirely directed upward ( Figure 2b) and reached about one third of the magnitude observed during the open water period. Due to increasing solar radiation, the air was heated rapidly in the second half of May, exceeding the freezing temperature as well as T 0 almost continuously since 21 May 2014. We observed an immediate switch to predominantly stable conditions and negative heat fluxes at this point, prevailing until ice-off. The start of the melting ice cover period was further indicated by a sharp R n increase around 21 May 2014, and the beginning of continuous changes in measurement height recognized by the snow depth sensor due to the melting of the ice beneath (but not below) the platform. During the short cold period from 27 to 29 May 2014, nearneutral conditions prevailed and the heat flux dynamics switched back to the dynamics observed during the frozen ice cover period, that is, net heat release. Snowfall on 29 May 2014 led to a period with a few very cloudy days (until 2 June 2014) with ongoing near-neutral conditions and very low turbulent heat fluxes, related to a low ΔT, that is, small differences between T a and T w0:5 m . Increasing T a caused a switch back to the conditions observed during the start of the melting ice cover period, that is, predominantly stable atmospheric conditions and downward H and LE. The latter indicated condensation or resublimation (deposition) at the lake (ice) surface. After the snowfall on 29 May 2014, which increased the albedo significantly, the albedo dropped to a level very similar to the albedo of open water, reflecting the crusting of the snow and the degradation of the ice.
As soon as the lake surface was ice free, LE became positive and continued as almost consistent evaporation with peaks of up to 274.8 W/m 2 , corresponding to 5.2 mm/day evaporative loss, and only a few short exceptions of condensation (Figure 2b). In comparison, H alternated between upward and downward fluxes depending on ΔT and atmospheric stability. Based on daily mean values, the mean and cumulative E for the open water period were 0.85 mm/day and 98.5 mm, respectively.
We observed only small differences in flux magnitudes and dynamics by comparison of "shore" (U 50°-150°) and "lake" (U 151°-49°) footprints ( Figure S2b; Table 1) during the period when the platform was close to the eastern shore of the lake. Footprint analysis indicated that during this period the lake contributed at least

10.1029/2017JD027751
Journal of Geophysical Research: Atmospheres 80% to the cumulative flux ( Figure S2c). Mean H and LE were higher than during the rest of the open water period, independent of the footprint.
The Bowen ratio β (Table 1) indicated that during the frozen ice cover and open water period, evaporation was a more important process of heat exchange at the surface than H. In comparison, during the melting ice cover period, median average deviation of β was high and H dominated the heat exchange.

Lake Thermal Dynamics
The thermal structure of the lake revealed a weak inverse stratification of the water column below the ice cover ( Figure 2f). The bottom of the water column remained unfrozen in winter with temperatures >2.7°C. We observed significant melting at the beginning of the melting ice cover period (22-26 May 2014) indicated by a pronounced vertical shift in the 0°C depth (about 35 cm). This shift was partly caused by the relative lift of the thermistor chain away from its initial position due to the melting of the ice cover surface beneath the platform and the associated relative elevation of the platform. Furthermore, the midday increase of T w surf to 0°C reflected the melt-freeze cycles on the ice cover surface, which stopped for some days starting with the 3-day cold period end of May. Water column temperatures already increased toward the end of the ice cover period and water column mixing, although not fully, had already started some days before complete ice-off. Soon after ice-off the thermal stratification broke down completely and the lake became isothermal. Since then, the lake warmed up rapidly reaching about 13°C on 11 July 2014, facilitated by the high absorption of solar radiation due to small albedo. The warming of the lake closely followed T a dynamics and was augmented by continuous and mostly full water column mixing sustained throughout major parts of the open water season. Stability indices (Lake Number and Wedderburn Number, not shown here) indicated that the water column stratified only four times for a few days during the open water period. During these short phases we observed low turbulent water-air heat exchange.

Energy Balance Components and Closure
The mean diurnal courses of the four EB components changed during the study period ( Figure 3). Magnitudes of R n and ΔQ were consistently higher than the turbulent heat fluxes. R n peaked around 2 pm (up to 320 W/m 2 during the melting ice cover period) and was negative during the night (down to

10.1029/2017JD027751
Journal of Geophysical Research: Atmospheres À41 W/m 2 during the frozen ice cover period). The largest daily amplitudes of R n and ΔQ were observed during open water conditions with 338 W/m 2 and 405 W/m 2 , respectively. The monthly averaged diurnal course of ΔQ showed a temporal shift of peak heat gain and loss from frozen ice cover toward open water conditions of about 4 and 6 hr, respectively. The diurnal course of ΔQ during the open water period was controlled by R n , with water column warming during the day, enhanced by thorough water column mixing. We expect the calculated ΔQ in the water column to be slightly overestimated by the shifting depths of the T w sensors due to the elevation of the raft and the adoption of T w0:5m as T w surf until the platform was floating.
The fluxes of H and LE had different mean diurnal patterns for the distinct lake surface conditions. During frozen ice cover conditions we observed very small turbulent heat fluxes during the night and consistently positive heat fluxes during daytime, peaking around 3 pm. In contrast, during the melting ice cover period, H and LE were consistently negative for the whole day, with the highest downward fluxes at around 5:30 pm. During open water conditions, the averaged fluxes were positive throughout the day with average daily maxima observed at around 6.30 am for H and 3.30 pm for LE. Heat lost from the water column was particularly transferred into H, indicated by the increasing H in the morning. Highest daily amplitudes of H and LE were observed during melting ice cover (56 W/m 2 ) and frozen ice cover conditions (44 W/m 2 ), respectively. The standard deviation (std) of the mean daily courses was the highest during daytime, and its dynamic most often followed the one of the mean diurnal course of the EB components (except, e.g., H during the open water period). The mean std for the average day of each subperiod was highest for ΔQ (Table S3).
The annual heat budget of the lake, as total amount of heat necessary to raise T w from minimum (0.87°C; 3 May 2014) to summer maximum (13.84°C; 1 August 2014), was 272 MJ/m 2 (see also Figure 2d), accounting for about 28% of the available radiant energy.
During the open water period the cumulative turbulent heat fluxes of an average day accounted for 98% of the available energy, corresponding to half-hourly averages of Res of about 1.5 W/m 2 . In comparison, during the ice cover period much more energy was available than spent, resulting in Res of 34.1 W/m 2 and 201.0 W/ m 2 for the frozen and melting ice cover period, respectively, accounting for 632 MJ/m 2 in total. Assuming that this Res corresponded entirely to the melting energy, we calculated an initial ice thickness of 2.06 m, compared to our observations of about 1.7-2.0 m at the beginning of the study period. The average ice decay rates during frozen and melting ice cover period were 0.96 cm/day and 5.67 cm/day, respectively.

Parameterization of Heat Fluxes
The controls UΔT and UΔe explained 58-91% of the heat flux variations during the three subperiods ( Figure S3). For both H and LE, the adjusted coefficients of determination (R 2 adj. ) for the linear relationships of the turbulent heat fluxes and their controls reflected best correlations for the open water period and particularly for H. In case of LE, excluding the 3-week period when the platform was close to the eastern lake shore (10-30 July 2014) from the data set led to slightly higher (R 2 adj. = 0.94) and lower (R 2 adj. = 0.67) correlations for H and LE for the open water period, respectively. Considering the whole study period UΔT and UΔe explained 81% and 73% of the variation in H and LE, respectively.
We parameterized the turbulent heat fluxes with one in situ and three independent models and evaluated their performance based on Taylor diagrams (Taylor, 2001;Figure 4; in addition Table 2 for the mean biases of modeled H and LE and Table S2 for R 2 adj. ). We found very similar results considering the whole study period and particularly the open water period during both stable and unstable conditions, for simplicity (z À d)/L ≥ 0 and (z À d)/L < 0, respectively. The model results were in moderate agreement with the measurements during the open water period, with modeled heat flux cycles that were sufficiently correlated with the observed (Pearson correlation coefficient of about 0.8). In contrast, the centered root mean square error (CRMSE, 20-30 W/m 2 ) and mean bias were nonnegligible considering the average measured heat fluxes during the open water period (Table 1). The high CRMSE was attributable to consistently underestimated heat flux amplitudes by all models, that is, std of the model results was much smaller compared to the measurements, indicating a strong underestimation of the small-scale temporal variability of the modeled fluxes. During the frozen ice cover period, moderate model results were achieved for both LE in stable and unstable conditions and H in unstable conditions, with a lower performance from IN-SITU v1 and SEA-ICE1. Differences between the models are mainly represented by std. In contrast, measured and modeled cycles of H were not correlated under stable conditions independent of the model. ΔT was exclusively positive (Figure 3b)

10.1029/2017JD027751
Journal of Geophysical Research: Atmospheres platform-related bias in T 0s (section 2.7) and therefore could not represent stable conditions. For the melting ice cover, instead, the strongest performance was found for H in stable conditions, with lower CRMSE and higher correlations for LAKE1 and SEA compared to IN-SITU v1. However, mean biases of the fluxes modeled were particularly high for H in these conditions, similarly to those during the frozen ice cover period. Whereas all models performed very poorly in case of LE during stable conditions, they showed no correlation to measured and modeled cycles of LE and negative correlations for H in unstable conditions. This seems to be the effect of an underestimated T 0 during daytime. The std of the modeled fluxes was most often lower than for the measurements, indicating an almost consistent underestimation of the heat flux variability (amplitude) by the models during the whole study period, which is also indicated by the much smoother diurnal courses of the model results compared to the measurements ( Figure S4). We found very small differences between model runs with fixed and dynamic measurement height indicating that the models are not sensitive to changes in measurement height of around 0.8 m. We yielded very similar parameterization results as the ones discussed for three additional parameterization schemes (LAKE2, LAKE3, and SEA-ICE2; Table S2).
We further evaluated the performance of IN-SITU v1, LAKE1, SEA, and SEA-ICE1 by calculating the EB for each subperiod with the modeled H and LE. During the open water period the EB was closed by 56%, 67%, and 59% and the back-calculated initial ice thickness was 186, 155, and 158 cm for IN-SITU v1, LAKE1, and SEA, respectively. This indicated an overestimation of H and LE by the models during the ice-cover period (particularly pronounced for H during the frozen ice cover period) and, in contrast, an underestimation during the open water period (particularly pronounced for LE), resulting in high mean biases ( Table 2). The daily courses of the modeled H and LE were very similar between the models and only slightly shifted (about 1 hr) compared to the measurements, except the more pronounced shift of about 4 hr in LE during the open water period ( Figure S4). The overall phasing (cycles) of the modeled heat flux dynamics was shown by Pearson correlation coefficient (Figure 4) to be comparable to the measured dynamics (exceptions mentioned above).
The effective transfer coefficients (C H and C E ; see Table S2, coefficients not available for SEA-ICE1 and LAKE3) of the independent models reflected their dependence on atmospheric stability. The models derived significantly higher coefficients during unstable than stable conditions, as heat loss is enhanced when the atmosphere is unstable. They are on average lowest during ice cover melt. For the schemes LAKE1 as well as LAKE2 and SEA C H and C E were congruent, respectively. In comparison, the in situ model (IN-SITU v1) revealed differences in C H and C E especially during ice cover conditions (C H > C E ), whereas the coefficients were very similar during the open water period. For the frozen ice cover period C H and C E were higher for SEA-ICE2 (1.83·10 À3 and 1.88·10 À3 ) than for the other models (mean C H, E = 1.59·10 À3 ).

Magnitude of Lake-Atmosphere Heat Exchange
We observed significant turbulent heat fluxes during the frozen ice cover period, which were, however, lower than during the open water period, as less incoming radiation and higher surface albedo limit the additional available energy and ΔT is smaller (Kirillin et al., 2012). The half-hourly heat flux measurements with clear diurnal cycle reflect the heat budgets of the thick ice and the (only temporary) snow cover rather than the water column underneath the ice, as ice and snow cover decouple the water column from atmospheric heat fluxes (Bengtsson, 2012). However, considering longer timescales, the course of the lake surface-atmosphere heat fluxes can be linked to the water column due to its interaction with the ice cover. Whereas the decoupling of the water column from the atmosphere by the ice cover maintains the persistence of its warmth (which is further maintained by the heat flux from the sediment to the water column, see, e.g., Bengtsson, 2012 andBoike et al., 2015), the warmer water column continuously affects the ice temperature by a small sensible heat flux to the underside of the ice (1-2 W/m 2 during winter in most lakes, Bengtsson, 2012; 5 W/m 2 in a boreal lake, Leppäranta et al., 2016). This small heat flux and the relatively high thermal conductivity of ice facilitates the persistence of heat fluxes to the atmosphere during winter on the long term (Jeffries et al., 1999). Furthermore, the lake was only rarely covered and thus insulated by a complete snow cover and prone to strong snow drift. The significant heat fluxes during the frozen ice cover period stand in contrast to very small heat fluxes observed above tundra during this time of the year, where the soil and small ponds freeze completely and the subsurface heat reservoir is emptied quickly at the beginning of winter (Langer et al., 2016).

Journal of Geophysical Research: Atmospheres
During the melting ice cover period we observed strong downward H and LE. At this time of the year, the tundra is most often already snow free and thus releasing heat to the colder atmosphere (Langer et al., 2011). Relatively warm and moist air is most likely advected over the still ice-covered lake, resulting in pronounced negative H and LE and initiating the heating of the surface and condensation or deposition. We expect condensation to dominate, as T 0 was about 0°C.
During the open water period, measured H reflected the atmospheric stability, with small fluxes during stable stratification due to suppressed turbulence. Potter (2011) reported for a thermokarst lake in Northern Alaska daily average H and LE for the open water period of 25.2 and 33.8 W/m 2 , compared to 14.1 and 54.9 W/m 2 during our study, respectively. The magnitudes of turbulent heat fluxes during the open water period were similar to those observed for the tundra (Langer et al., 2011;Ohmura, 1982), contradicting Jeffries and Richter-Menge (2015), reporting that seasonally ice-covered lakes have the highest evaporation rates of any high-latitude terrestrial surface during ice-free periods. According to Rouse et al. (2005), strong solar radiation absorption and water column heating, as well as the continuously wet surface result in high lake evaporation rates during the open water season. The evaporative loss during our study period was about one fifth lower than the average annual sum of rainfall (about 125 mm; see Boike et al., 2013) and is at the lower limit of the estimated annual evaporative losses of 0.4 to 2.3 mm/day for thermokarst lakes in the Arctic Coastal Plain during early summer (Arp et al., 2015;see also Potter, 2011). Arp et al. (2015) observed that annual evaporation rates are mainly controlled by the lake ice regime, with bedfast ice lakes showing higher evaporation rates than floating ice lakes due to earlier ice-free conditions. We missed about 1.5 months (late summer/early autumn) of the evaporative season in our campaign, which is potentially characterized by large convective heat fluxes due to the high cumulative heat storage and the cooling overlying air (Petrone & Rouse, 2000;Rouse et al., 2005).

Lake Heat Storage and Water Column Mixing
Starting already at the end of the ice cover period, the lake was an efficient heat absorber and distributed heat effectively, first through convective and, following ice-off, additional wind-driven mixing. The warming of the water column due to solar radiation penetrating the ice cover, further indicated by diurnal changes in T w below the ice, was already visible during the frozen ice cover period in the absence of a consistent snow cover. According to Bengtsson (2012), penetration of solar radiation through the snow-free, thinning ice cover in spring can account for 10-40 W/m 2 (see also Leppäranta et al., 2016). Absorption is probably enhanced by increasing ice surface wetness and below-ice circulation and convective mixing of the water column, which can be controlled by the heat flux from the sediment toward the end of the ice cover melt (Bengtsson, 2012;Kirillin et al., 2012;Jammet et al., 2015). Ice breakup and ice-off in 2014 were within the average timing and the initial ice thickness was similar as reported for other thermokarst lakes in the central Lena River Delta (Boike et al., 2015).
Starting at ice-off, it took about 18 days to warm the water body to an almost consistent level, with only slight changes up to when we stopped taking measurements. The almost continuous mixing of the water column is consistent with Boike et al. (2015), stating that only arctic lakes with mean depths >5 m show continuous summer stratification of 1 month or longer. The annual heat budget of the lake is in the range reported by Boike et al. (2015) and Rouse et al. (2005) for lakes of similar depth and size.

Energy Balance Closure and Melting Energy
The open water EBC (98%) supports the reliability of the EC measurements and supports our initial assumption that EB components other than R n , ΔQ, and turbulent heat fluxes played a negligible role. For lakes with mean depths of 1-5 m Boike et al. (2015) modeled a maximum heat flux of about 3-4 W/m 2 into the sediment during the summer and about 0-7 W/m 2 from the sediment into the water column during the ice cover period (see also Bengtsson, 2012). We expect the heat flux due to precipitation to be negligible, as precipitation is low and dominated by light rain events in the central Lena Delta . Net heat flux through lake inlets and outlets as well as groundwater is assumed to be negligible as, for example, the contribution of surface runoff and suprapermafrost groundwater is generally much smaller compared to vertical water balance components (precipitation and evaporation). The EBC is not commonly estimated and discussed in lake studies, especially including EC measurements. The EBC calculated for the open water period of the studied 10.1029/2017JD027751 Journal of Geophysical Research: Atmospheres lake is at the upper limit of EBC reported for a subtropical lake and the open water period of two boreal lakes of 84-98% (Liu et al., 2009), 72-82% (Nordbo et al., 2011), and 70-99% (Mammarella et al., 2015), respectively.
The diurnal cycle of ΔQ with changing phase shift during the study period, supports our decision to calculate Res and EBC for an average day of each month. However, the seasonal heat storage change is not considered. We expect strong release of the heat, which was added to the lake within our study period, during late summer and autumn before refreeze (see Petrone & Rouse, 2000;Rouse et al., 2005), a period we did not cover with our measurements. High Res for the ice cover period (frozen: 34.1 W/m 2 ; melting: 201.0 W/m 2 ) confirmed that the energy required for ice cover melt represents an important EB component (Lund et al., 2017). The back-calculation resulted in a slightly overestimated initial ice thickness but an ice decay rate during the melting ice cover period which is comparable to other studies (e.g., 1.5 cm/day for boreal lake, Jakkila et al., 2009; 6-10 cm/day for floating ice lakes, 11-16 cm/day for bedfast ice lakes in Northern Alaska, Arp et al., 2015). Other potential reasons for the energy imbalance during the ice cover period include the decoupling of processes of the water column and the atmosphere (Kirillin et al., 2012), biases in the simplified correction of the net outgoing radiation, an overestimation of ΔQ s due to the platform elevation, the different footprint of the sensors measuring the EB components (McJannet et al., 2011;Vercauteren et al., 2009;Wilson et al., 2002), and the nonconsideration of snow cover in the balance.

Evaluation of Heat Flux Parameterization Schemes
As land surface models implemented in numerical weather prediction and climate models are sensitive to transfer coefficients, the adoption of in situ derived coefficients as reported in this study can improve the estimation of the surface heat fluxes (Deng et al., 2013;Potter, 2011;Xiao et al., 2013). During the open water season the in situ and LAKE1 model-derived C H10 and C E10 are within the range of values estimated based on the open water season measurements at other lakes of different depths, sizes, and climatic regions (e.g., Mammarella et al., 2015;Xiao et al., 2013). Potter (2011) derived mass transfer coefficients C H10 and C E10 for the ice-free period of a thermokarst lake in Northern Alaska of 2.303 W·m À3 ·s·°C À1 and 35.293 W·m À3 ·s·kPa À1 , respectively. Accordingly, forcing the intercept of the regression to be 0, we calculated C H10 and C E10 of 1.591 W·m À3 ·s·°C À1 and 42.446 W·m À3 ·s·kPa À1 for the Siberian Lake, respectively.
An important source of error in bulk flux algorithms is the accuracy of the input variables (Bourassa et al., 2013). The in situ parameterization and model comparison revealed that the model results are less sensitive to the choice of the specific model than to the accuracy of the surface water temperature measurement (see also Brunke et al., 2002;Fairall et al., 1996). For example, we showed that an overestimated T 0s can result in a reverse temperature gradient, failing to represent downward H in the model results for stable conditions during the frozen ice cover period (particularly during the night; see Figure 3). All models performed very poorly for H and LE during unstable conditions when the ice cover was melting, which seemed to be caused by underestimated T 0m during daytime and the associated bias (reversion in case of H) of UΔT and UΔe.
Modeled H during the open water period followed similar, slightly shifted diurnal dynamics compared to our measurements, but on lower magnitudes especially during daytime, and switching to negative values (as also indicated by negative ΔT for about 6 hr in the afternoon; Figure 3f). The reason for this discrepancy is not as obvious as for the frozen ice cover period, as we expect measurements of T 0m to be accurate during this period. We further exclude flow distortion and an associated higher proportion of upward winds and overestimated H to be the reason, as we filtered the fluxes appropriately (Text S1). Woolway et al. (2015) have discussed the risk of underestimating turbulent fluxes across the air-water interface for lakes by classical bulk transfer models for oceans, as they do not account for lake-specific conditions such as wind variation caused by sheltering and fetch limitation of wind (for example, Schladow et al., 2002). Furthermore, the atmospheric boundary layer stratification over lakes is more variable due to stronger heating and cooling by the adjacent land surface and lower U (Verburg & Antenucci, 2010). The strongly underestimated diurnal variability and magnitude of LE as well as the remarkable shift in the diurnal cycle compared to the measured fluxes during the open water period, which cannot be explained by data inaccuracies, raise concerns on the ability of the models to simulate lake evaporation.
The independent parameterization schemes developed for open water conditions performed similarly well during the frozen ice cover conditions as the sea ice-specific scheme. Sea ice-specific models consider long-life stably stratified atmospheric conditions, for example, which often develop above snow and ice 10.1029/2017JD027751 Journal of Geophysical Research: Atmospheres surfaces (Andreas, 2002;Andreas et al., 2010;Atlaskin & Vihma, 2012) and distinct physical properties of the ice pack and the surface including roughness length (Brunke et al., 2006;Lu et al., 2013;Persson et al., 2002). However, during the frozen ice cover period, where we tested the refined models, we observed most often near-neutral conditions.
Apart from the sensibility of the models to water temperature accuracy, the performance of the models was mostly constrained by the underestimation of the heat flux variability. However, in any case, care should be taken when comparing measurements of EC fluxes against bulk flux estimations. The former represent single point measurements, while the latter are ecosystem-scale integrated estimations with some assumptions on the vertical stratification and its effect on the vertical turbulent transport. The single-point measurements are affected by turbulence intermittency, and the time averaging of point measurements does not always provide identical values to the ensemble averaging of turbulent fields. Besides, such factors as varying footprint conditioned by the occasional movement of the measurement platform, inevitably affect the EC flux measurements. Therefore, the EC results cannot be considered as an ultimate benchmark for bulk models, although they can provide a useful insight into their ability to simulate surface fluxes.

Summary and Conclusions
We have measured and modeled H and LE of an arctic Siberian thermokarst lake with the aim of contributing to the understanding of their dynamics during three basic lake states: frozen ice cover, melting ice cover including ice breakup, and open water. The lake showed strongly differing heat flux dynamics for the three subperiods in terms of magnitude, direction, and temporal patterns. The consistent significant H and LE release during the frozen ice cover period highlighted the importance of lakes for the landscape heat budget and water balance during winter. Negative H and LE during melting ice cover conditions indicated local heat advection and condensation on the melting ice cover surface. Along with enhancing incoming radiation, these downward heat fluxes provided the energy for ice cover melting, which was highlighted as a very important component of the energy balance during spring. In contrast, during the open water period an EB including net radiation, lake heat storage change, and the turbulent heat fluxes was sufficiently closed.
The measured heat fluxes reflected the atmospheric stability conditions and followed the bulk transfer relationships during the three lake states, supporting the reliability of our EC measurements. However, we found inaccuracies in the measured and calculated surface water temperatures that originated from our measurement setup. These inaccuracies constrained the ability to parameterize the turbulent heat fluxes based on the bulk transfer relationships during the ice cover period. In addition, strongly underestimated diurnal variability and magnitude as well as the shift in the diurnal cycle of modeled LE compared to the measurements raise concerns on the ability of the models to simulate evaporation of an arctic thermokarst lake during the open water period. For all three subperiods, an in situ and three independent models performed similarly and differences were most often attributable to the reflection of the heat flux variability. We could not support the necessity of applying model algorithms that are specifically adapted to ice cover conditions. The in situ derived transfer coefficients can contribute to a better parameterization of the physical processes of arctic lakes and the land surface interaction in the Arctic. The applicability of these coefficients and the independent models needs, however, to be further investigated. Future field campaigns should focus on the ice cover period and consider the practical challenges in a well-conceived measurement setup, to improve the accuracy of the measurements and with it the performance of in situ parameterizations and the evaluation of independent models.