2.1 Study site

The field campaign was carried out on 21–27 May 2014 in Kilpisjärvi (Fig. 1) – a “mid-size” (surface area: 37.1 km²) lake located in the Scandinavian Mountains in northwestern Finland (69^∘03^′ N 20^∘50^′ E, 473 m above sea level). The average and maximum depths of the lake are 19.5 and 57 m, respectively. The average ice-on and ice-off dates are 9 November and 18 June (mean value over 1952–2015), the earliest recorded ice-on date was 21 October and the earliest and latest breakup dates were 2 June and 1 July, respectively (data of the station of the Finnish Environment Institute (Korhonen, 2006). In summer, surface water temperatures reach 10–15 ^∘C in mid-August. On average, the maximum annual ice thickness is reached in mid-April, amounting to 89 cm (the range over 1952–2015 was 77–114 cm). Weather data (standard observations, including air temperature at 2 m height and wind speed at 10 m height) were adopted from the Kilpisjärvi Biological Station and Enontekiö Kilpisjärvi Kyläkeskus (EKK). The latter station was founded in Kilpisjärvi in 1962 by the Finnish Meteorological Institute. The annual mean air temperature is between −4 and −1 ^∘C, with mean monthly air temperatures below 0 ^∘C from October to April. Snow thickness on the ground in March–April is on average 90 cm. The snow thickness on the lake ice tends to be lower than on the surrounding land, partly due to snow removal by winds, and partly due to transformation of snow to snow ice (Leppäranta et al., 2018), with an annual maximum of 37±9.1 cm.

Figure 1(a) Geographical location of the study site, (b) bathymetric map with the position of the thermistor chain, and (c) the grid used in the model of free oscillations. Blue arrows show the river flow direction.

Download

2.2 Vertical heat transport in the water column

The background information on the vertical temperature distribution was collected by a thermistor chain installed close to the lake center (${\mathrm{69}}^{\circ }{\mathrm{2}}^{\prime }{\mathrm{3.09}}^{\prime \prime }$ N, ${\mathrm{20}}^{\circ }{\mathrm{46}}^{\prime }{\mathrm{45.87}}^{\prime \prime }$ E, depth 32 m). The chain was equipped with nine TR-1060 temperature loggers (RBR, Canada, declared accuracy 0.002 ^∘C, sampling period 10.0 s) at heights 31.25 27.00 25.00, 20.80, 18.65, 14.40, 8.25, 6.20, and 4.15 m above the bottom, and a RBR-DUO temperature–pressure logger (RBR, Canada, temperature accuracy 0.002 ^∘C, pressure accuracy 0.05 % of 50 dbar, sampling period 2.0 s) at 1.00 m above the bottom. Another temperature–pressure logger recorded atmospheric pressure at the lake shore for correction of the water-level records. Several casts of vertical profiles of temperature and conductivity were performed during the study using the RINKO CTD profiler (JFE Advantech, Japan, temperature accuracy 0.01 ^∘C, conductivity accuracy 0.01 mS cm⁻¹, sampling rate 10 Hz). Ice and snow thickness were measured during the CTD casts with an accuracy of ±0.5 cm.

Data on photosynthetically active radiation (PAR) were collected by three compact PAR loggers DEFI (Alec Electronics and JFE Advantech, Japan), one on the ice surface, one frozen into the ice cover, and one in the water column at 0.3 m beneath the ice cover. All sensors recorded planar quantum irradiance with 10 min intervals. The spectrum of solar short-wave (wavelength range 200–2500 nm) radiation is strongly modified by lake waters: clear water (or ice) quickly absorbs the long-wave (infrared) part of the spectrum and yellow substance absorbs the short-wave (ultraviolet) part. As a result, at <1 m depth, > 95 % of the penetrated radiation falls within the PAR spectral range of 400–700 nm (see e.g., Jerlov, 1976). In humic brown-water lakes, like Kilpisjärvi, the equivalence between PAR and total solar radiation is even closer (e.g., Leppäranta et al., 2010). Hence, we adopted the measured PAR values as characteristic of the total downward short-wave radiation flux. Transformation between the measured quantum irradiance I_q (µmol s⁻¹ m⁻²) and the net downward short-wave radiation Q_I (W m⁻²) was obtained in the atmosphere by the relation $I_{\mathrm{q}}/Q_I=\mathrm{4.60}$ µmol J⁻¹ (see Leppäranta et al., 2010). In water and ice the coefficient depends on the color and turbidity of water. In Finnish and Estonian lakes the transformation coefficient has been found to vary between 4.8 and 5.5 µmol J⁻¹, larger values referring to more turbid waters (Reinart et al., 1998). For the clear-water Lake Kilpisjärvi the coefficient was chosen as 4.8 µmol J⁻¹. The accuracy of the PAR sensors is ±2 W m⁻² in air and ±3 W m⁻² in water or ice. For a comprehensive description of PAR measurements, see Leppäranta et al. (2018). Light extinction within the water column was estimated from the vertical profiles of PAR taken during the previous field campaign in 2013 (Kirillin et al., 2015). The profiles were measured in 2013 with spherical quantum PAR sensor LI-193SA (LiCor, USA) attached to a CTD-90M profiler (Sea & Sun Technology, Germany). The light extinction coefficient γ and the radiation value at the ice–water interface Q_I0 were determined from a one-band exponential approximation of the short-wave radiation profile Q_I(z) in the water column,

The vertical turbulent heat flux within the bulk of the water column Q(z,t) as a function of time t and depth z was estimated from temperatures measured by the thermistor chain T(z,t) using the “flux-gradient method” which adopts the one-dimensional equation of heat transfer, with horizontal advection neglected:

where $C_{\mathrm{p}}\mathit{\rho }\approx \mathrm{4.18}\times {\mathrm{10}}^{\mathrm{6}}$ J K⁻¹ m⁻³ is the product of the water heat capacity and density. The solar radiation flux Q_I(z,t) was recovered from PAR measurements and Eq. (1). Integration of Eq. (2) from a reference depth H, usually chosen close to the lake bottom, to a depth z yields the turbulent flux Q(z,t) as

Assuming that, at the end of the ice-covered period, the bottom sediment is in thermal equilibrium with the deep water layers, i.e., the vertical turbulent flux Q(H)≈0, Eq. (3) was solved numerically using finite differences for differentiation and the trapezoid method for integration.

2.3 Measurements in the ice–water boundary layer

Measurements of the turbulent pulsations of current velocities were performed using pulse-coherent high-resolution acoustic Doppler profiler HR-Aquadopp (Nortek AS, Norway). The profiler was attached to a thin plate with positive buoyancy and deployed alongside the ice cover base, with the acoustic head oriented looking downwards at a right angle to the instrument. The profiler registered three components of the current velocity vector at a sampling rate of 0.5 Hz with a vertical resolution of 0.02 m.

The HR-Aquadopp data were used for estimation of the dissipation rate ε of the turbulent kinetic energy (TKE). The procedure of ε estimation generally followed that from the previous study (McGinnis et al., 2015): the velocity structure function

was calculated from the HR-Aquadopp data and was used for estimation of the TKE dissipation rate, as described by Wiles et al. (2006). Here, v(z) is the along-beam velocity fluctuation at distance z, r is the depth range of the dissipation estimation, and angle brackets denote time averaging. Three estimations of the velocity fluctuations were determined by extracting the mean velocity value for the averaging periods of 10, 20, and 30 min. Three values of the maximum estimation range for the velocity correlation r=0.4, 0.5, and 0.6 m were tested, covering the range recommended by Wiles et al. (2006) for weakly stratified turbulence. The TKE dissipation rate was estimated by fitting the equation

Here, the constant C_v is 3^1∕3 (see, e.g., Lien and D'Asaro, 2002). The fitting constant Noise representing the average effect of the acoustic noise on the velocity fluctuations was used for the goodness-of-fit check, using the condition

The measurements, for which Eq. (6) was valid, were abandoned as noisy. Further quality check was performed by comparison of the 27 arrays of the TKE dissipation rate calculated from the three HR-Aquadopp beams with three different values of r and three averaging periods. The discrepancy between the different estimations did not exceed 10 % (not shown). The estimations based on the averaging time of 20 min were adopted for the further analysis averaged over the three beam estimations. The lower value of r=0.4 m was adopted to avoid occasional influence on the results of the increased instrumental noise at higher distances from the acoustic head.

To measure the vertical distribution and short-term variability of water temperature within the strongly stratified ice–water boundary layer, a custom measurement setup was used. It consisted of a temperature registration platform with 12 T-Solo fast16 fast-response temperature loggers (RBR, Canada, accuracy 0.002 ^∘C, time constant 0.1 s, sampling frequency 16 Hz) equipped with stainless steel sledges on the upper side and attached to remote operated vehicle VideoRay Pro4 (VideoRay, USA). The distance from the ice base to the uppermost sensor was 5 cm; other loggers were distributed vertically at distances of 2 cm between the sensors to a depth of 27 cm from the ice base. The platform was able to slide along the ice bottom up 150 m away from the deployment site to locations undisturbed by instrument installation. Two loggers of the downward solar radiation were additionally attached to the platform, but malfunctioned during the deployment. Series of 24 to 48 h long deployments were performed with continuous temperature registration at a sampling frequency of 2 Hz.

2.4 Two-dimensional model of standing barotropic waves

The variations in the temperature and turbulent mixing within the ice–water boundary layer were analyzed for their relationship with the standing barotropic waves (seiches) by using spectral analysis and modeling of the free oscillations of the lake. Spectral densities were determined using the standard periodogram method. The free oscillation modes of Kilpisärvi were estimated from numerical solution of the eigenvalue problem

on the geometry of the lake in the horizontal coordinates x and y. Here, $c=(gH{)}^{\mathrm{1}/\mathrm{2}}$ is the long-wave speed (celerity) based on flat-bottom approximation of the lake bathymetry with the mean depth H, λk is the eigenfrequency and Ψ is the eigenfunction of the harmonic level oscillations ζ in time t, defined as

Equation (7) was solved on a triangular mesh approximating the Kilpisärvi geometry (Fig. 1c) using the finite element solver of COMSOL Multiphysics software. The lake geometry was provided by National Land Survey of Finland.

2.5 Approximate solutions of the vertical heat transport equation

As mentioned above, the ice–water interface in small lakes is commonly assumed to be turbulence-free. Assuming the heat transport within IL is steady-state ($\partial Q\partial t^{-\mathrm{1}}=\mathrm{0}$) and non-turbulent ($Q\left(z\right)=-\mathit{\kappa }\partial T\partial z^{-\mathrm{1}}$), a direct solution of the steady-state heat transport Eq. (2) turns to

Subject to boundary conditions (Fig. 2),

Equation (9) yields the solution (Barnes and Hobbie, 1960)

Here, δ is the IL thickness, κ is the heat conduction coefficient, $I=(C_{\mathrm{p}}\mathit{\rho }{)}^{-\mathrm{1}}{Q}_I$ is the kinematic flux of the short-wave solar radiation, and T_m is the temperature of the convectively mixed layer beneath IL. The IL thickness δ can in turn be found by applying the condition of smooth temperature transition between IL and the thermally homogeneous CL beneath, i.e.,

or

The latter equation can be solved with respect to δ if the function I(z) is known. For a single-band exponential decay of the short-wave solar radiation within the water column (1), the solution (12) turns to

This equation is solved numerically with respect to δ. As it follows from Eq. (13), δ depends on the CL temperature T_m, the radiation flux I₀ and the water transparency γ: the more turbid the lake is, the stronger the radiation is, the higher T_m is and the thinner IL is. The water–ice heat flux is determined directly from the solution (10) as

i.e., the upward heat flux equals the solar radiation absorbed within the IL in the steady-state non-turbulent conditions and no temperature gradient (zero heat flux) at the IL bottom.

Another simple model, particularly useful for the analysis below, is that of non-penetrative growth of the CL on the background of stable density stratification in the quiescent layer (QL) beneath. The approximation was first introduced by Zubov (1945) and was later named “encroachment” (Carson and Smith, 1975) in contrast to the “entrainment” regime, which involves entrainment of a denser fluid from QL to CL. If the well-mixed convective layer develops on top of the QL by increase in the mixed layer temperature T_m, without changing the quasi-linear temperature gradient Γ in the QL, the rate dh dt⁻¹ of the CL deepening and the heating rate are related as

which suggests zero turbulent (heat) buoyancy flux at the CL bottom. Together with the previous assumption of the zero flux at the CL–IL boundary, the CL development can be estimated from data on solar radiation flux:

This regime is strictly never realized in natural convection (Carson and Smith, 1975), but can be closely approached, at early stages of convection development (Carson, 1973).

Figure 2Schematic of the temperature structure of the under-ice water column illustrating the major assumptions under models (10)–(14) and (15)–(16).

Download

3.1 Weather and ice conditions

During the observation period of 22–27 May, the air temperature did not change much, with the mean value of 2.15 ^∘C ±1.75, minimum of −2.00 ^∘C and maximum of 5.30 ^∘C. Low winds of 2–4 m s⁻¹ prevailed most of the period, with a moderate wind increase up to 9 m s⁻¹ on 22 May and a short event of strong winds > 12 m s⁻¹ on 24–25 May (Fig. 3). Conditioned by mostly cloudy weather and a thick ice cover with more than 50 % of snow ice, the short-wave radiation flux right beneath the ice cover amounted to daily means of 10–20 W m⁻², with daily peak values of 35–70 W m⁻² (Fig. 3), which made up 6 %–7 % of the solar radiation flux at the ice surface.

Figure 3Wind (bars) and under-ice radiation (solid line) during the observation period.

Download

The ice was snow-free (snow thickness < 1 cm). The ice thickness at the beginning of the field campaign was close to the seasonal maximum values: mean 97 cm, min 82 cm, max 126 cm, based on 29 measurements from a cross section along a short axis of the lake; 52 cm of the ice cover consisted of “white ice” (snow ice) or frozen slush on top of the congelation ice. From 20 to 26 May the ice thickness reduced to 73 cm (minimum 61 cm, maximum 80 cm, based on 37 measurements from a cross section along a short axis of the lake).

For the period under investigation, the ice surface temperature remained close to the freezing point of 0 ^∘C, which also implies that the entire ice cover had a nearly constant temperature of 0 ^∘C. The heat budget of the ice cover is dominated in this case by the heat exchange at the ice bottom and by the absorption of solar radiation within the ice cover, whereas the conductive heat transport across the ice cover is small. The situation is characteristic of the early stage of the ice cover melt and was considered in detail by Aslamov et al. (2014). Since half of the ice sheet was snow ice (white ice), its low transmittance limited the transfer of sunlight to water and the ice also experienced internal deterioration from absorption of solar radiation. At the upper boundary the heat balance was positive and caused surface melting. The surface heat balance was dominated by the solar and long-wave radiation balance, turbulent fluxes being small. Upper surface melting and internal melting accounted for 1–2 cm per day. Further details on the structure, optical properties, and heat budget of the ice cover are presented elsewhere (Leppäranta et al., 2018).

3.2 Temperature structure of the water column

The vertical temperature distribution had a three-layer structure, typical of the spring convection in ice-covered lakes due to solar heating (Fig. 4; cf. also Fig. 2): the well-mixed convective layer (CL) had initial thickness h_m≈5.5 m and temperature T_m=0.8 ^∘C, adjacent to a 1.5 m thick stratified interfacial layer (IL) on top with temperature decrease to the freezing point at the ice–water interface. The quiescent layer (QL) beneath the CL was nearly linearly stratified with downward temperature increase in Γ≈0.04 K m⁻¹ (hereinafter Kelvin is used for temperature gradients, while absolute temperatures are given in degrees Celsius for convenience). Solar radiation penetrating the ice cover caused an increase in temperature in the CL at $\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{2}}$ ^∘C day⁻¹ and an increase in its thickness at 0.33 m day⁻¹.

Figure 4Temperature profiles at the beginning (solid line) and at the end (dotted line) of the field study. The inset shows a zoomed portion with CL and temperature gradients at its boundaries.

Download

In contrast to observations at later stages of the convection development in ice-covered lakes (see, e.g., Mironov et al., 2002), there was no appreciable temperature/buoyancy jump at the bottom of the CL: the well-mixed layer h developed on top of the QL by increase in the mixed layer temperature T_m, without changing the quasi-linear temperature gradient Γ in the QL, as described by the “encroachment” model (Eq. 15). At the top of the CL, in turn, convection produced an increase in the temperature gradient across the IL from 0.58 K m⁻¹ on 20 May to 1.74 K m⁻¹ on 26 May that corresponds to an increase in the conductive heat flux across the IL from 0.33 to 0.99 W m⁻² according to the model of Barnes and Hobbie (Eq. 14).

3.3 Heat budget within the water column

The light attenuation within the water column was estimated based on 24 PAR profiles taken in May 2013 as $\mathit{\gamma }=\mathrm{0.30}\pm \mathrm{0.05}$ m⁻¹. We adopted this estimation for the 2014 campaign, in suggestion of low year-to-year variability of γ for the same season. The resulting vertical profiles of the radiative flux (Fig. 5c, d) were used for estimation of the vertical turbulent heat transport (black lines in Fig. 5a, b) by vertical integration of the heat equation (Eq. 3). During daytime, the turbulent fluxes were generally directed downwards, increased across the convective mixed layer towards the lake surface to ∼10 W m⁻², and dropped to almost zero within the stratified surface layer at the ice base, an indication of the turbulence damping by stratification in the surface layer. Exact values of the turbulent fluxes in the vicinity of the ice–water interface were difficult to estimate by the heat budget method due to coarse vertical resolution. No appreciable negative fluxes existed at the bottom of the convective mixed layer (∼7 m depth), in agreement with the absence of an entrainment layer in the temperature profiles. Also, the low but non-zero turbulent fluxes below the convective mixed layer (at 10–25 m water depths) may have resulted from the method constraints, but could be produced by natural processes, such as horizontal heat advection of vertical mixing by internal waves. During nighttime, turbulent transport changed its direction to upward. The fluxes were generally below 10 W m⁻², except the night of 24–25 May, when a strong flux of ∼18 W m⁻² developed near the ice base, decreasing downwards across the mixed layer.

Figure 5Vertical turbulent (a, b) and radiation (c, d) heat fluxes in the bulk of the water column. (a, c) show the daytime averages; (b, d) correspond to the nighttime values.

Download

3.4 Dynamics of the stratified interfacial layer (IL)

The stratified IL occupied the up to 2 m thick layer under the ice base, separating the CL from the ice base. High-resolution temperature observations during 24–26 May in the IL demonstrated periodic fluctuations with the 24 h period and an amplitude of 0.2–0.3 K, related to the subsurface heating by the penetrating solar radiation. In the background of the diurnal cycle, higher-frequency temperature oscillations persisted during the first half of the observation period (Fig. 6a). Irregular high-amplitude temperature fluctuations lasted for the first 5–6 h (24 May, 13:00–19:00 in Fig. 6a), and were replaced later by harmonic quasi-sinusoidal oscillations (24 May 19:00–25 May 19:00 in Fig. 6a) decaying within ≈24 h. Concurrently with the strong temperature fluctuations, the TKE dissipation rate ε increased by an order of magnitude $>{\mathrm{10}}^{-\mathrm{8}}$ W kg⁻¹ in the entire IL, as estimated from the velocity fluctuations (Fig. 6b). The burst in the turbulent mixing decayed within 5–6 h to the background values at the limit of the detectable turbulence level of $\mathit{\epsilon }\approx {\mathrm{10}}^{-\mathrm{10}}$–10⁻⁹ W kg⁻¹. The burst in TKE dissipation and temperature fluctuations coincided with the wind increase on 24–25 May (Fig. 3).

Figure 6(a) Temperature fluctuations in the IL and (b) TKE dissipation rates in the upper water column. In white areas of (b) no developed turbulence could be detected, because no inertial subrange could be identified according to Eq. (5).

Download

The dataset on the TKE dissipation rate covered a longer period than the temperatures in the IL, and demonstrated a close connection between wind speeds over the lake and the turbulence intensity under the ice cover: apart from the strong wind burst on 24–25 May, the winds of ∼10 m s⁻¹ on 22 May increased ε to $\sim {\mathrm{10}}^{-\mathrm{8}}$ W kg⁻¹ (cf. solid line and wind bars in Fig. 7). The whole series of 22–27 May suggests that wind has an appreciable effect on mixing under ice at speeds > 10 m s⁻¹; otherwise, the mixing intensity in the IL remains low, with $\mathit{\epsilon }\approx {\mathrm{10}}^{-\mathrm{9}}$ W kg⁻¹. Weaker increases in ε up to 10^−8.5 W kg⁻¹ took place at quasi-diurnal intervals and were connected apparently to the periodical increase in radiation-driven convection in the CL underneath (cf. temperature variations in Fig. 6a).

Figure 7Vertically averaged TKE dissipation rate under ice (solid line) vs. wind speed over the lake (bars) and vertical temperature gradient from the two uppermost temperature loggers under ice, at 5 and 7 cm under the ice bottom (line with symbols). Single squares show the mean daily vertical temperature gradient under ice according to the steady-state model of Barnes and Hobbie (1960); see Eq. (14). Note the different vertical axes for temperature, wind, and ε.

Download

The vertical temperature gradients $\partial T\partial z^{-\mathrm{1}}$ under ice (line with circles in Fig. 7) demonstrated strong correlation with dynamics of wind and ε during the same period, following the decay of the mixing intensity from its peak on 24 May. Vertical conductive heat flux, calculated from these gradients as $-\mathit{\rho }{c}_P\mathit{\varkappa }\partial T\partial z^{-\mathrm{1}}$, remained as low as 1.0–1.5 W m⁻². However, the high dissipation rates under ice (Fig. 6b) suggest the heat transport within this layer remains mainly turbulent and exceeds the conductive one. The estimations of the water–ice heat flux from the model of Barnes and Hobbie (1960, see Eq. 14) are approximately 2 times higher (2.0–3.0 W m⁻²), but do not reflect the temporal variability connected to ε, because they rely on the assumption of a steady-state IL.

Spectral analysis of the under-ice temperature and pressure records demonstrated a strongly periodic character of water motions. The spectral energy of water-level (pressure) variations (Fig. 8a) was concentrated at several peak frequencies, fairly corresponding to the eigenfrequencies from the two-dimensional seiche model (Eq. 7). However, the three different stages distinguished in the subsurface temperature time series (cf. Fig. 6a) – “wind burst” on 24 May, “free oscillations” on 25 May, and “convection only” on 26 May – reveal rather different spectral patterns: during the “wind burst”, most of the energy still resided at the low-frequency “gravest seiche mode” oscillations (triangles in Fig. 8a). The energy peaks in temperature oscillations during the same period (Fig. 8b) were relatively low, suggesting that the vertical motions in the stratified IL were only weakly periodic. The “free oscillations” stage was characterized by concentrations of motions around the frequencies of the second and higher seiche modes, i.e., vertical motions with the maximum amplitudes away from the lake shores. During the “pure convection” stage the spectral energy in the range of short-term oscillations 1–100 cph was essentially zero (not shown); i.e., at relatively low winds of <5 m s⁻¹, no seiche oscillations were produced under the ice cover.

Figure 8Power spectral densities for (a) water-level variations and (b) temperatures in the IL for the “wind burst” on 24 May 2014 12:00–24:00 (red triangles) and “free oscillations” on 25 May 2014 03:00–15:00 (black circles). The vertical dotted lines mark the modeled frequencies of free oscillations in (flat-bottom) Kilpisjärvi.

Download