Modelling circulation in an ice-covered lake

In deep ice-covered lakes with temperatures below 4 °C the heat flux from the bottom sediment results in a horizontal density gradient and a consequent flow along the bottom slope. Measurements in Lake Pääjärvi, Finland, show a stable temperature field where a heat gain through the bottom and a heat loss through the ice nearly balance each other. The circulation is thermal with low velocities (less than 1.5 cm s). We used the 3D hydrodynamic Princeton Ocean Model as a tool to simulate the water circulation and the temperature distribution under the ice. The model forcing was based on field temperature measurements. The model simulations suggest that in midwinter the velocity field of the upper water layers is anticyclonic while that of deep layers is cyclonic. Comparison with current measurements at one site showed good agreement between the modelled and observed results. On the basis of the modelled results it is possible to better understand the distributions of some microorganisms and the accumulation of oxygen depleted waters in the deepest part of the lake.


INTRODUCTION
Thermal structure and circulation dynamics of lakes have a major impact on lake ecosystems.Due to seasonal and synoptic weather variations, the research of this impact is a challenging task.In the boreal zone this is especially true since in winter ice cover radically reduces the forcing by atmospheric factors (Adams 1981;Leppäranta 2009).Although some basic features of winter thermo-and hydrodynamics are known, our understanding is predominantly based on theoretical aspects and heuristic reasoning rather than field measurements (see Salonen et al. 2009).In the 1940s Mortimer found that winter warming of water in deep lakes has an important role in under-ice hydrodynamics (Mortimer 1941;Mortimer & Mackereth 1958).Bengtsson (1996) presented four different hydrodynamic processes, which can cause currents under ice: river through-flow, wind-driven elastic waves in the ice cover, penetration of solar radiation and sediment heat release.The heat flux from the sediment increases the density of overlying cold water and causes slow benthic water currents along the bottom slope.
Present technology enables us to collect precise data with a high spatial resolution.This, combined with modelling, allows us to go beyond the classical static approach towards the real description of dynamic processes in ice-covered lakes.With models we can describe and quantify the water flow and temperature fields and examine the effect of different forcing.In addition, with well-calibrated models we can reliably obtain information without measurements on-site.
In model application the determination of boundary conditions is important and we have to focus on the most important fluxes.The in-and outflows are important in lakes with short retention times.For a lake with snowcovered ice the solar radiation forcing is small and when the water exchange is slow, the boundary heating determines the circulation.In this study these fluxes are handled by fixing the boundary temperature.For a more advanced treatment, ice as well as sediment sub-models would be needed.
In this study we hypothesize that in medium-sized ice-covered lakes density gradients induced by bottom heating generate a lake-wide circulation system.Our main tool is a model application for Lake Pääjärvi in Southern Finland based on the 3D hydrodynamic Princeton Ocean Model.We also use field data to determine the boundary temperature field for the model and to validate the results.First we apply the hydrodynamic model to an idealized conical basin and then extend the model to Lake Pääjärvi.The results show the importance of Coriolis acceleration and the high stability in the winter-time circulation of medium-sized lakes.

Site
Lake Pääjärvi is a medium-sized, deep, dimictic, humic and mesotrophic lake (Fig. 1).It is located in the boreal vegetation zone (61°04′N, 25°08′E), its surface area is 13.4 km 2 and its mean and maximum depths are 15 and 87 m, respectively.The average retention time is 3.3 years, but during winter due to low river discharges the water exchange is much slower than average.The water is brown in colour due to a high humus content and the light attenuation coefficient is 1.5 m -1 (Arst et al. 2008).
The long-term ice climatology of Pääjärvi described by Kärkäs (2000) showed no significant changes during the 1900s.On average, the ice season starts within the period 25-30th of November and ends within the period 20-25th of April.During the 2003-2004 winter the ice cover formed on 21 December and melted by 28 April.The maximum annual ice thickness is on average 45 cm, while in 2004 it was 41 cm.In winter months there is normally 10-15 cm of snow on top of the ice.

Measurements
We collected water temperature and current velocity data in winter 2004.Water temperature was measured between the littoral zone and the deepest part of the lake (Fig. 1).Temperature measurements were made along two transects on 12 February 2004 with a CTD (FSI MCTD-3, Falmouth Scientific, Inc.), which according to the manufacturer has an accuracy of ± 0.005 °C for temperature, ± 0.005 mS cm -1 for conductivity and ± 0.02% for pressure.The western line comprised of 17 and the middle line of 13 vertical profiles with a distance of 50-200 m between the measurement points.Water temperature near the deepest site of the lake was measured four times during midwinter.
In winter 2004 continuous water current measurements (Shirasawa et al. 2006) were conducted with an acoustic small-scale flow meter in the western bay of the lake at a depth of 5 m (Fig. 1).A three-dimensional electromagnetic current meter (ACM 32M, Alec Electronics Co., Ltd., Japan) was used to record the currents at 30-min time intervals.The velocities were low, not much above the measurement threshold, and the results can be taken to be reliable for their magnitude, directionality and stability.The measurement period began before ice formation and extended beyond ice break-up in April, thus providing a good illustration of the role of the ice cover in the winter circulation.

Modelling
The Princeton Ocean Model (Blumberg & Mellor 1987) was used to estimate the under-ice circulation of Pääjärvi.This model is a nonlinear, fully three-dimensional, primitive equation, finite-difference model including the Mellor & Yamada (1982) level 2.5 turbulence submodel.For horizontal diffusion the Smagorinsky diffusivity is used.The model is widely used for ocean and lake circulation research and also for studying boundary layer problems.Two basins were examined.Firstly, a simplified conical lake with a hypsographic curve corresponding to Pääjärvi was modelled to find the general circulation features in such basins.Secondly, Lake Pääjärvi bathymetry was used.Schwab et al. (1995) had used a simplified conical lake bathymetry to manifest the cyclonic circulation in deep stratified lakes during the open water period.
In the model the -coordinate σ system was used with 30 vertical -levels.

σ
We had a maximum layer thickness of 2.8 m in the deepest part of the lake and a layer thickness of 0.5 m on average.The horizontal size of grid cells was 50 x y ∆ = ∆ = m.In total the Lake Pääjärvi grid had 179 × 107 cells in the horizontal direction.The basin with an ideal circular geometry had a grid with 94 × 94 cells with a horizontal size of 50 In the present case the flow regime is basically nonturbulent, low Reynolds number flow.It is known that in such conditions the Mellor-Yamada turbulence model results differ from the observed values (Baumert et al. 2000).However, with the limited validation data available this model seemed to produce realistic effective heat conductivity (see Eq. 1) and drag force so that the general momentum and heat transfer characteristics as well as the overall circulation level were realistic.
The forcing in both basins was taken care by fixed boundary conditions for temperature.This means that the bottom and surface temperatures were constant during the simulation but temperatures in the other layers of the water body were solved by the model.The initial temperature field was based on data collected on 12 February 2004 from transects over the lake and assuming 0 °C temperature just beneath the ice cover (Fig. 2).The later CTD sections showed that in February-March the temperature distribution was quite stable, and therefore the results represent the diagnostic circulation field under a snow-covered ice lid over the lake.A model run was finished when the system reached equilibrium with the existing boundary conditions.

Water temperature and heat fluxes
Temperature measurements in the lake show winter stratification with an inverse thermocline at about 5 m depth and bottom temperatures of 2-3 °C (Fig. 2).The upper 5 m layer was slightly warmer in the centre than in the peripheral part of the lake with isotherms sloping down towards the shores, but horizontal differences were small in deeper water (Fig. 2).Near the bottom a slight opposite sloping of isotherms was observed along both transects, which was more pronounced along the western transect.
The stratification was stable through the winter with a heat flux from the bottom compensating for the heat loss through the ice.At the deepest site the water temperature profiles revealed a consistent vertical thermal structure over the midwinter 2004 observation period.Water temperature profiles showed persistent thermoclines in 10-30 m and 55-75 m layers.At the bottom of the deepest site the water temperature was up to 0.5 °C higher than in the water column above.
The approximation of heat fluxes through the bottom BL ( ) and surface boundaries SL ( ) is difficult in icecovered lakes since the water flow is in the laminar or laminar-turbulent transition regime.Turbulent flows may occur occasionally in mid-winter and more so after snow melt due to forcing by solar radiation.The boundary heat fluxes can be estimated with a semiempirical approach as (Petrov et al. 2006; see also Jakkila et al. 2009 where w,eff K is called the effective heat conductivity.When the flow approaches the laminar regime, w,eff w , K K → which is the molecular conductivity of heat.The molecular conductivity is 0.6 W/(m °C) but the effective conductivity is one order of magnitude larger.
Here it is assumed that the molecular approximation can be used at the bottom layer (Table 1).In the surface layer the same is true when the snow cover is thick enough to prevent the transmission of solar radiation through the ice.This is probably an underestimate for the level of heat and momentum transfer but serves as a good first approximation.
The obtained heat flux from the bottom to the water was smaller than the flux from the water to the atmosphere through the ice.However, both fluxes were quite small so that the net loss was not significant enough to be reflected in the temperature profiles.In the model simulations below we assume that under-ice circulation is controlled by the bottom heat flux.

Current measurements
The observations show that water currents at 5 m depth had slow velocities (mean velocity 0.85 cm s -1 ) in the western bay (Fig. 3).The time series shows a highly stable structure for the current speed and direction, but some variations were observed.In January the currents were turning from the north towards the east.In February the decreasing currents turned back towards the north and remained at a low velocity level during March.This may indicate a local decreasing heat flux during the winter.The mean current over the fourmonth period was directed towards the deep middle part of the lake (Fig. 3).
According to the vertical water temperature distribution, the thickness of the bottom boundary layer was approximately 3 m.If we take the mean water temperature in this layer to be 2.75 °C and the observed mean velocity of 0.85 cm s -1 , the Reynolds number is 1.6 × 10 4 .Since the current speed ranged by 50% around the mean, we can see that the flow was just within the transition zone towards the turbulent regime.

Simulation of winter circulation in a conical basin
The initial water temperature field of the idealized conical basin was chosen with higher temperatures at the slopes than in the central part corresponding to observations in Lake Pääjärvi (Fig. 4a).The resulting isotherms are convex across the lake.The initial current velocities were set to zero.The model simulation conserved the convexity of the isotherms in other depths except at 30-40 m due to adjustments towards geostrophic balance (Fig. 4b).
In the under-ice lake circulation, inertial, advective and frictional forces are small, and therefore the Coriolis acceleration becomes an important factor in balancing the pressure gradient.Indeed, due to the thermal structure of the water body, an anti-cyclonic gyre was dominating the flow field in the upper layers (Fig. 5a, b), while in deeper layers the circulation was cyclonic (Fig. 5c-f).

Simulation of the winter circulation in Lake Pääjärvi
In the Lake Pääjärvi simulation the initial temperature field was similar to that in the conical basin showing higher temperatures at the slopes (Fig. 6).The basic features in the simulation results were as in the conical basin, with the anticyclonic and cyclonic gyres resulting from the balance between the pressure gradient and the Coriolis acceleration.Comparison of the mean current velocity (Fig. 3) and the model result for the steady state circulation (Fig. 7a) shows agreement between the direction and magnitude (measured mean velocity at 5 m depth vs. model velocity at 2 m depth).In both cases the direction of the current was towards the northeast with a speed of 0.5-1.0cm s -1 .It should be noted that    the observed current was not in a steady state but was slowly varying.As discussed above, the boundary drag forces are likely to be biased downwards.However, this is not a major issue since friction is of minor importance in the present situation.In the simulation the initially convex isotherms were modified at 30-40 m depths as in the conical basin and likely for the same reasons.
Using the molecular heat conductivity at the boundaries also underestimates the heat fluxes but this underestimation holds both for the gain at the bottom and for the loss at the surface, and therefore the vertical temperature profile is consistent with the observations as is shown in Fig. 8.But it is clear that the interaction of the water body with the boundaries is biased in the model.
The modelling results demonstrate that to reach the steady state from no motion for fixed boundary conditions, a one month model run is needed.The results are shown in Fig. 7a-f.In this case an extensive cyclonic gyre with a maximum velocity of ca 10 mm s -1 was formed in the surface layer, covering the area from the southwestern shore towards the northwestern shore and continuing to the northern shore and then back to the southern shore (Fig. 7a).At 10 m depth two gyres were formed in the centre of the lake.They were converging on the northeastern side of the deep basin, with a current towards the deepest site of the lake, where they merged with the waters within the cyclonic gyre (Fig. 7b).A somewhat stronger and spatially extensive cyclonic gyre appeared at 20 m depth (Fig. 7c).The same pattern was found in deeper levels with diminishing velocity (Fig. 7d, e).

DISCUSSION AND CONCLUSIONS
Temperature structure and circulation in an ice-covered lake were investigated considering field data and results from a 3D hydrodynamic model.Field work was made in Lake Pääjärvi, southern Finland, while modelling concerned this lake and an idealized conical basin corresponding in dimensions to Pääjärvi.
Both for the conical basin (Fig. 5) and Lake Pääjärvi (Fig. 7) the simulation produced two overlying gyres with a changing flow direction at intermediate depths.This suggests that in medium-sized ice-covered deep lake basins with a similar bottom heat storage capacity, anti-cyclonic circulation in the upper layer and cyclonic circulation in the lower layer are formed.This pattern is forced by the density gradient induced by sediment heat release and influenced by the Earth's rotation or the Coriolis acceleration.The necessary conditions for the formation of these gyres are the heat storage capacity of the lake bottom sediment and lake morphology.It is not clear how large a deviation from the conical basin geometry is required to change the circulation characteristics qualitatively.
Sinking of water along the slopes to the deeper layers bends the isotherms to follow the bottom topography.Upwards directed flow of water at the lake centre is due to this double-gyre structure.The circulation pattern is similar to the atmospheric tropical cyclone or to Ekman pumping induced by the atmospheric cyclone (Gill 1982), although the nature and magnitude of forcing are totally different.
Many large lakes and semi-enclosed marginal seas and estuaries show persistent net cyclonic circulation (e.g.Cummins & Foreman 1998).This pattern is most apparent during the period when the lakes are stratified.Velocities in this circulation are generally low compared to the transient water currents induced by wind.However, because of its persistent character, cyclonic circulation can contribute significantly to the thermal structure in the lake and the redistribution of dissolved and suspended materials in the water body.In Pääjärvi the effect of the Earth's rotation during summer stratification has been demonstrated earlier (Virta 2001) but its role is much less important than in winter.The winter conditions are quite different as the ciculation is thermal, while in the open water season it is predominatly wind-driven.
Several explanations have been proposed for the predominantly cyclonic circulation in open stratified basins.Emery & Csanady (1973) suggested a mean cyclonic curl in the wind stress field as a possible mechanism.Schwab et al. (1995) considered, as in this study, that this circulation in large stratified lakes is a result of isotherm shape near the lake bottom, causing internal pressure gradients resulting in geostrophic cyclonic circulation.
We found that the Princeton Ocean Model is a suitable tool for the prediction of under-ice water circulation as it contains terrain following -coordinates, σ and well developed turbulence closure.It has been observed in other lakes that the fluxes change in time and in space (e.g.Thanderz 1973).However, in our application the initial temperature field and also the near-bottom water temperature used during the simulation were determined on the basis of only a few measurements.Using fixed boundary conditions is somewhat unrealistic, since the sediment releases heat to the water and it cools down.A likely explanation of the goodness of the Princeton Ocean Model in the present application is that the circulation is governed by the pressure gradient and the Coriolis acceleration for which the model provides a very good representation although the -coordinates σ may cause a small erroneous component on geostrophic flows due to numerical truncation error in the calculation of the pressure gradient (Mellor et al. 1994(Mellor et al. , 1998)).We consider this error small and used the -coordinates σ approach as it allows the simulation of the bottom boundary layer, which was necessary for this study.
A recent study by Jakkila (2007) (see also Jakkila et al. 2009) provides some indications of temporal and spatial variations in sediment heat fluxes in Pääjärvi.Jakkila measured the heat exchange between the water and atmosphere in detail.Using the temperature profiles measured in different sites of the lake, he estimated that the change in heat storage of the lake increased by ~ 2 W m -2 in February-March 2004.The maximum increase in the heat storage (6 W m -2 ) was observed in the deep middle parts of the lake, which indicates an advection of heat from the littoral zone to the deepest parts of the lake.
Our observations and simulations clearly show that in ice-covered lakes cooling below 4 °C causes a sediment heat flux, which generates a significant circulation regime.Although the flow velocities are rather slow, their persistency over the long winter makes them important factors for micro-organisms.Also, the displacement of sediment oxygen consumption towards the deepest water layers may have far-reaching biological consequences.
In the future we shall concentrate on the dynamic simulation of under-ice currents including solar and atmospheric forcing and also extend the simulations to the period of ice-out.We will also include dynamic values of the sediment heat fluxes.

Fig. 1 .
Fig. 1.The location (a) and bathymetric map (b) of Lake Pääjärvi.The surface area of the lake is 13.4 km 2 , volume 0.2 × 10 9 m 3 and maximum depth 85 m.Black dot -the location of the current meter; circle -the deepest site; solid line -transect in model output; dashed lines -transverse cross sections in the western bay (left) and middle basin (right).

Fig. 2 .
Fig. 2. Water temperature isotherms (°C) for transverse cross sections of the lake across the western bay (upper panel) and the middle of Pääjärvi (lower panel) on 12 February 2004.

Fig. 3 .
Fig. 3. Water currents at the western bay station of Pääjärvi between 1 January and 1 April 2004.(a) Direction and velocity distribution, (b) time series of velocity magnitude and (c) direction of the current clockwise from the north.

Fig. 4 .
Fig. 4. The temperature fields in the simplified conical lake case: (a) the initial temperature distribution (°C, in mid-February 2004), (b) the simulated water temperature (°C) in steady state (y-axis 80 m).

Fig. 5 .
Fig. 5.The horizontal velocity distribution at different depths in the simplified conical lake case.

Fig. 6 .
Fig. 6.The temperature fields in the real Pääjärvi simulations: (a) the initial water temperature distribution (°C, in mid-February 2004), (b) the simulated water temperature (°C ) in steady state in the Pääjärvi model transect (y-axis 80 m).

Fig. 7 .
Fig. 7.The horizontal velocity distribution at 2-40 m depths (a-e), and depth-averaged velocity distribution (f) in the real Pääjärvi case.The location of the current meter is indicated with a black dot and the deepest site with a circle.

Fig. 8 .
Fig. 8.The observed (°C, in mid-February 2004) and the simulated water temperature (°C) in steady state in the Pääjärvi model at the deepest point of Lake Pääjärvi.

Table 1 .
Estimation of bottom and surface heat fluxes based on temperature profiles at the deepest site of Pääjärvi