Introduction

Observations of freeze-up date, maximum ice thickness, and break-up date of lakes in middle- to high-latitude regions have been suggested as useful indicators of changes in local and regional climate (Reference Palecki and Barry.Palecki and Barry, 1986: Reference SchindlerSchindler and others, 1990; Reference Robertson, Ragotzkie and Magnuson.Robertson and others, 1992). Since frozen lakes provide an index of integrated seasonal temperature trends, the analysis of the seasonal evolution of frozen lakes provides a valuable contribution to climate monitoring, particularly in the data-sparse and climate-sensitive high-elevation and high-latitude regions of the world.

In this paper an energy-balance model describing the seasonal evolution of lake ice is developed and used to address the sensitivity of lake-ice growth and decay mechanisms under variations in atmospheric forcing. Model output is validated against lake-ice observations made during the winter of 1991/92 from Great Slave Lake in northern Ganada and during the winter of 1992/93 in Glacier National Park, northern Montana, U.S.A. This model has been found to simulate successfully the growth and decay features of lakes in Glacier National Park (Reference Liston and Hall.Liston and Hall, 1995).

The one-dimensional, unsteady model is composed of four major sub-models. First, a lake-mixing, energy-transport sub-model describes the evolution of lake-water temperature and stratification. Ice is initiated when the upper lake temperature falls below freezing. Secondly, a snow sub-model describes the depth and density of snow as it accumulates, metamorphoses, and melts on top of the lake ice. Thirdly, a lake-ice growth sub-model produces ice by two different mechanisms: relatively bubble-free clear ice grows at the ice/water interface due to thermal gradients within the ice, and snow ice forms at the lake-ice surface from the freezing of water-saturated snow, or slush: this slush can form from the upwelling of water due to an overburden of snow, from snowmelt or from rain on snow. Fourthly, a surface energy-balance sub-model is implemented to determine the surface temperature and energy available for freezing or melting. A schematic illustrating the components of the model and a representative temperature profile is found in Figure 1. In addition to complete energy-balance components over the annual cycle, key model output includes the dates of lake freeze-up and break-up and the end-of-season clear-ice, snow-ice, and total ice depths. At a minimum, the model is driven with daily atmospheric forcing of precipitation, wind speed and air temperature. It is also capable of taking advantage of observed humidity and radiation components.

Fig. 1.

Schematic illustration of key features of the lake-ice growth model and a representative temperature profile.

General Model Equations

A detailed description of the model equations and the solution procedures can be found in Reference Liston and Hall.Liston and Hall (1995). A more general survey of the model will be presented in this section.

Lake model

To determine the initiation of ice formation on the lake, a model is required which describes the temperature evolution of the lake during the late summer, autumn and early winter cooling periods. To accomplish this, consider the following one-dimensional heat-transfer equation,

(1)

where T _W is the water temperature, t is time, z is the vertical coordinate measured downward from the water surface, K _H is the turbulent diffusiviiy for heat, k _W is the water thermal conductivity, ρ _W is the water density, and C ^P is the specific heat of water. The radiation flux, q, is given by

(2)

where Q ₀ is the net solar radiation penetrating the surface, and η is the extinction coefficient, assumed to be 0.6 m⁻¹ (Reference AshtonAshton, 1986).

The vertical distribution of the difrusivity is strongly a Function of surface wind shear stress and local stratification. Following Reference Henderson-Sellers.Henderson-Sellers (1985), the diffusivity is given by

(3)

where the Richardson number for the water body, Ri_W, is

(4)

the surface shear velocity, w _s*, is

(5)

and the shear velocity decay constant, k*, is given by

(6)

where k is von Karman’s constant, U ₁₀ is the wind speed at 10 m, ϕ is latitude, and P ₀ is the turbulent Prandtl number at neutral stability, assumed to be unity.

The Brunt–Vaisala frequency, N ², is determined by the density stratification, which requires a description of the water density variation with temperature (Reference Kell. and Franks.Kell, 1972). If a computed temperature profile is unstably stratified, the model produces an instantaneous mixing which eliminates the unstable configuration.

To solve the lake-transport equation, initial conditions and surface and lower boundary conditions must be provided. For an initial temperature condition, the water column is assumed to be isothermal at 4°C. In the absence of ice, the surface temperature is determined from the surface energy balance; if ice is present, the temperature of the water/ice interface is equal to the water freezing temperature. The heat flux from the lake bottom is recognized to play an important role in governing the thermal regime of small, shallow lakes. For the relatively large and deep lakes considered in this paper, a zero diffusive-flux boundary condition is applied at the lower boundary. To initiate ice on the lake surface, a 1 mm thick ice layer is formed when the temperatures of the lake surface layers are less than the water freezing temperature.

Ice and snow models

A model describing the growth of clear lake ice can be developed by performing an energy balance at the ice/water interface. While the water/ice boundary is fixed at the stable-equilibrium freezing temperature, the freezing process produces latent heat which must be conducted through the ice. These processes, together with the convection occurring within the water, lead to a thermal energy balance at the ice/water interface that takes the form

(7)

where ρ _i is the ice density; L _f is the latent heat of fusion; h _W is the convective transfer coefficient; T _s0 is the surface temperature; T _W is the water freezing temperature; T _W is the lake water temperature; z is the depth of each layer of thermal conductivity, k, where the subscripts i, s, m and w indicate individual layers of ice, snow, a mix of snow and water, and water, respectively; and dz _i/dt is the velocity of the moving ice boundary. Since the lakes considered in this study are typically snow-covered, the influence of solar radiation penetrating the ice has heen assumed to be of negligible importance in this formulation.

Solution of this equation requires knowledge of snowpack depth and thermal conductivity. The effective thermal conductivity of snow, k _s, is defined to be a function of snow density, ρ _s, (Reference VerseghyVerseghy, 1991),

(8)

Density changes occur in the snowpack by two mechanisms. First, density increases can result from compaction, and follow Reference AndersonAnderson {1976),

(9)

where T _s is the snow temperature, W _s is the weight of snow expressed in water-equivalent depth, and A ₁ and A ₂ are empirical constants. The snow is modeled as a single layer, and the snow temperature is determined from the arithmetic mean of the surface and snow ice interface temperatures, where the temperature at the snow/ice interface is determined from energy-balance considerations under steady-state conditions. The second density-modifying process results from snowmelting. The melted snow reduces the snow depth and is redistributed through the snowpack until a maximum snow density, assumed to the 550 kg m⁻³, is reached. Any additional meltwater accumulates at the base of the snowpack as a layer of water-saturated snow, or slush.

The energy flux available for melt, Q _m (see below), is first used to melt any available snow and then the lake ice. Under freezing conditions, the freezing of water held within any slush mixture of snow and water adds to the lake-ice depth in the form of snow ice.

Given input of liquid-equivalent precipitation, the precipitation is assumed to fall as snow if the wet-bulb temperature is lower than 1°C (Reference RogersRogers, 1979) . Precipitation falling as rain contributes to a liquid-water store on the lake surface. Snow falling on a bare ice surface or an unsaturated snow surface accumulates as new snow. Following Reference AndersonAnderson (1976), the new snow density, ρ _us, is given by

(10)

This new snow is added to the existing snowpack and the snow depth and mean density values are updated.

A further mechanism by which water is added to the snow cover is from the upwelling of lake water through cracks in the ice. This is the result of the inability of the ice to support the overburden of snow cover at a level where the top of the ice equals that of the water, A floating body displaces a volume of liquid equivalent to its weight. Applying this balance to the system of snow and water lying above an ice cover which is depressed to a point where the ice top meets the lake surface, leads to the buoyancy balance

(11)

where each term represents the force per unit area, or pressure, when multiplied by the gravitational acceleration, g. When the downward force (righthand side) exceeds the upward force (lefthand side), cracks in the lake ice are assumed to form or be present which allow lake water to saturate the ice surface. Under this condition, the snow cover adds to the upward buoyancy force, and the water that had previously accumulated on the ice surface is interconnected with the lake water. This leads to a new force balance

(12)

which can be solved for the new storage of water on the lake surface.

(13)

This water store is available to form snow ice depending on the available freezing energy.

Surface energy balance

The surface energy-balance equation is

(14)

where Q _si is the solar radiation reaching the surface of the Earth. Q _li is the incoming long-wave radiation, Q _le is the emitted long-wave radiation, Q _h is the turbulent exchange of sensible heat, Q _e is the turbulent exchange of latent heat, Q _c is the energy transport due to conduction, Q _m is the energy flux available for melt, and α _s is the albedo of the surface. Complete details of the equations comprising each term in this energy balance are provided by Reference Liston and Hall.Liston and Hall (1995).

The solar radiation term takes into account the time-dependent solar irradiance at the top of the atmosphere, the solar elevation angle and the net sky transmissivity. The downward long-wave radiation under clear skies and standard atmospheric conditions is given by Reference SatterlundSatterlund (1979). In the current study no attempt has been made to modify this formulation for Q _li to account for the presence of clouds. The long-wave radiation emitted by the various surfaces is computed under the assumption that they emit as grey bodies. The turbulent exchange of sensible and latent heat follow Reference Price and Dunne.Price and Dunne (1976). The conductive heat flux is given by the product of the thermal conductivity and the temperature gradient at the surface.

The roughness length varies with surface type, and in this application is assumed to be 0.001, 0.005 and 0.000001 m over open water, snow and bare ice, respectively (e.g. Reference EaglesonEagleson, 1970). The albedo, α _s decreases with increasing snow density following Reference AndersonAnderson (1976):

(15)

for snow densities ranging between 50 and 450 kg m ⁻³. For densities greater than 450 kg m⁻³ the albedo is assumed to vary linearily according to

(16)

Lake ice is assumed to have an albedo of 0.25, and a take surface without snow or ice has an albedo of 0.06 (e.g. Reference EaglesonEagleson, 1970).

When coupled to the coupled lake, lake-ice and snow-models, the solution of the energy balance provides the surface temperature boundary condition which forces the seasonal evolution of lake-water temperatures, lake-ice growth and decay, and snow-cover accumulation and metamorphism. With the support of daily wind speed, precipitation, cloud-cover fraction and maximum and minimum air-temperature observations, the energy-balance equations are solved iteratively for the surface temperature. T _s0. using the Newton Raphsou method. In the presence of snow and/or ice, surface temperatures T _s0 > 0°C resulting from the surface energy balance indicate that energy is available for melting, Q _m. The amount available is then computed by setting the surface temperature to 0°C and recomputing the surface energy balance. A similar procedure is adopted to compute the energy available to freeze. Q _i, liquid water within the snowpack.

Model Simulations

The coupled lake-ice atmosphere model described above is validated by considering two lakes located in significantly different climatic regions: St Mary Lake in eastern Glacier National Park (GNP), Montana, and Great Slave Lake in Northwest Territories, Canada. Once the model has been shown to perform satisfactorily, it will be used to analyze the sensitivity of lake-ice features to variations in atmospheric forcing. Field observations of snow and lake-ice depth were obtained for St Mary Lake during the winter of 1992/93, and daily meteorological observations were collected at the St Mary GXP ranger station, located near the mouth of St Mary Lake. These meteorological data include daily maximum and minimum temperature, water-equivalent precipitation, cloud cover and wind speed, and cover the period September 1992–May 1993. Lake-ice and snow-cover data for Great Slave Lake during the winter of 1991/92, and air-temperature, precipitation, wind-speed, station-pressure, cloud-fraction, and relative humidity data, for the period July 1991–July 1992, were obtained from the Atmospheric Environment Service, Canada. The ice- and snow-depth measurements were made on Back Bay, about 100 m off of the old Wardair float base, near Yellowknife, Northwest Territories, and the meteorological data were collected in Yellowknife. On Back Bay, Great Slave Lake, the 1991/92 snow depths and lake-ice depths are both within 15 cm of the mean 1961 90 values. During the period 1958-90, maximum and minimum lake-ice depths vary by as much as 50 cm from the mean depth value (1961–90), and as such the period 1991/92 is considered fairly representative of the average ice and temperature conditions for Yellowknife/Back Bay. The periods of model integral ion coincide with the acquired atmospheric forcing data. For the purposes of the following discussions, the term “freeze-up” refers to the earliest date on which the water body was completely covered with ice, and “break-up” refers to the date on which the water body was finally clear of all ice.

Snow accumulation on a relatively flat surface, such as that of a lake, is dependent upon complex interrelationships between wind speed, air temperature and precipitation, in combination with historical factors which serve to define the shear strength of the existing snowpack. To account for the obvious connection between reduced snow accumulation in response to higher wind speeds, the precipitation accumulating on the lake is determined by multiplying the observed snow water-equivalent precipitation by the empirical scaling factor β, where β has been chosen to equal 0.5 for St Mary Lake, and 0.75 for Great Slave Lake (see the meteorological data discussion below). Analysis of the St Mary meteorological data shows that during October and November the average daily air temperature was below freezing on several occasions, and stayed below freezing through December and almost all of January. In late January and early February, air temperatures rose above freezing before dropping during the last half of February. March began with above-freezing temperatures and, after a below-freezing period in the middle of the month, the air temperature rose and slaved above freezing. The December–April air temperature, wind speed and precipitation averaged −3.1°C, 6.6 m s⁻¹, and l.l mm d⁻¹, respectively. Minimum temperature, maximum wind speed and precipitation during the same period were −34°C, 22 m s⁻¹ and 15 mm d⁻¹, respectively.

The Great Slave Lake meteorological record shows that the air temperature dropped below freezing in early October and rose above freezing in early May. The December–April air temperature, wind speed and precipitation averaged −19.3’C, 3.9 ms⁻¹ and 0.6 mm d⁻¹, respectively. Minimum temperature, maximum wind speed and precipitation during the same period were −43°C, 10 m s⁻¹ and 9 mm d⁻¹, respectively.

The ice and snow depths simulated by the model for St Mary Lake are compared with observations in Figure 2a. Also included are the mean, maximum and minimum observed depths for each observation time. The number of observations varies from one Held excursion to the next, and ranges from one to nine observations per trip to the lake. Measurements were not always made in the same part of the lake during each trip. The model simulation and the average of the observations agree quite well, and consistently lie within the range of the observed values. The freeze-up date is simulated within the limits of the observations and occurs during the last week in November. Frequently, when an observation indicates little or no variability, only one or very few measurements were made. The model produces a maximum total ice depth of 63 cm in late February. The relatively warm periods in early February and March have led to the reductions in ice depth shown. The final melting of lake ice begins in the middle of March and continues until the ice is entirely melted by 1 April.

The ice and snow depths simulated bY the model for Great Slave Lake are compared with weekly observations in Figure 2b. The model simulation indicates that ice began forming on Great Slave Lake (Back Bay) on 12 October and total break-up was achieved by 7 June, Observations show that ice began forming on 16 October and the lake was ice-free on 2 June. The maximum observed ice depth on Great Slave was 137 cm; this is also well simulated by the model. The maximum observed snow depth on the lake was 40cm, and that produced by the model was 32 cm. Both the shape of the snow- and ice-growth curves and their magnitude have been captured in the model simulations.

Fig. 2.

Comparision between the ice and snow depths simulated by the model and the observations for: a. St Mary Lake. b. Great Slave Lake.

To explore the sensitivity of lake ice to variations in climate, several changes in atmospheric forcing were considered for the two climatic regimes currently found in northern Canada and Montana. These include modifications to air temperature, ±4°C; precipitation (accumulation), ±20%; wind speed, ± 3 m s⁻¹; and cloud-cover traction, ± 3/10. These modifications were applied to the meteorological observations used in the St Mary and Great Slave Lake simulations described above.

Figure 3a illustrates the changes in ice depth for St Mary Lake that result from air temperature changes of ± 4°C. This lake is located in a climatic region for which air temperatures are often near freezing. As a result, there exists a strong lake-ice sensitivity to changes in air temperature, and, for the case of increasing temperatures, this led to an ice-free date that is 4 weeks earlier than the control simulation. Changes of ± 20% in accumulation of precipitation on the ice surface have modified the ice depths slightly and have had no appreciable effect on the timing of break-up (not shown). Variations in wind speed can significantly alter the surface temperature regime. Reductions (increases) in the wind speed have lowered (raised) the surface temperature and led to increased (decreased) ice depths (Fig 3b). The cloud-cover fraction affects the incident solar radiation, and influences the energy available to melt the snow and ice, particularly during the spring months when temperatures are near freezing and incoming solar radiation is becoming more dominant. For St Mary Lake, changes in cloud tener of ± 3/10 produce mid-March ice depth differences of 10 cm, or roughly 25% of the ice cover, and shift the ice-free date by up to 10 d (Fig, 3c).

Fig. 3.

Ice depth variations for St Mary Lake in response to changes in: a. air temperature, b. wind speed, c. cloud-cover fraction.

Modification of the air temperature by ±4°C for Great Slave Lake leads to changes of only 5 cm in the maximum ice depth (Fig. 4). For the case of a 4°C temperature decrease (increase), the initiation of melting has been delayed (advanced) by nearly 3 (2) weeks, and the ice-free date has been delayed (advanced) by 2 weeks. The reduced-temperature scenario has produced a 2 week advance in freeze-up date. Changes in accumulation have not affected the break-up date, but have produced differences in maximum ice depth of roughly 10 cm (not shown). Reduced wind speeds led to an increase in the end-of-winter ice depth of nearly 10 cm, and produced a surface temperature low enough to delay the onset of break-up by two weeks (Fig 4b). Increased wind speeds produced smaller differences in ice depth and break-up date. Cloud-cover increases (decreases) had virtually no impact on the ice depths, and led to a break-up date that was 1 week later (earlier) than the control (Fig 4c).

Fig. 4.

Ice depth variations for Great Slave Lake in response to changes in: a. air temperature, b. wind speed, c. cloud-cover fraction.

Summary And Conclusions

To assess the response of lake ice-related parameters to potential changes in atmospheric forcing, a physically based computational model of the coupled lake, lake-ice, snow and atmosphere system has been developed (Reference Liston and Hall.Liston and Hall, 1995). The model is driven, at a minimum, with atmospheric forcing of air temperature, wind speed and precipitation. In addition to providing complete energy-balance components over the annual cycle, model output includes the seasonal evolution of lake-water temperatures, lake freeze-up and break-up dates and the complete growth and decay cycle of clear ice and snow ice within the lake-ice matrix. Model performance has been validated for two different climatic regions using meteorological and lake-ice observations Great Slave Lake in northern Ganada (1991/92) and from St Mary Lake in Glacier National Park, Montana (1992/93). To explore the sensitivity of lake ice to variations in climate, several changes in atmospheric forcing were considered for the two climatic regimes currently found in northern Canada and Montana.

Model integrations suggest that air-temperature changes of ±4°C can delay or speed up the freeze-up dates and break-up dales by as much as 4 weeks for St Mary Lake, and 2 weeks for Great Slave Lake. For both lakes, break-up date is more sensitive to air temperature changes than is freeze-up date. Snow accumulation changes of ±20% produce end-of-winter ice-depth changes of approximately 5–10 cm for both lakes, with no change in the break-up date. The imposed changes in wind speed, ± 3 m s⁻¹, modify the maximum ice depth of the lakes by 5–10 cm. For Great Slave Lake, the reduction of wind speed led to a surface temperature low enough to delay the onset of break-up by 2 weeks. Changes of ±3/10 cloud-cover fraction indicate a shifting of break-up dates by roughly 1 week.

For climate changes occurring at the magnitudes considered in this study, the above simulations provide important information concerning the requirements of in situ and remote-sensing observations and instruments designed to detect climate changes using lake ice-related parameters. To detect the changes addressed in the above simulations, instrumentation to measure ice depth would require sufficient resolution to detect end-of-winter lake-ice depth differences on the order of a few centimeters, in order to detect the 5–10 cm changes in ice depth. Detection of the changes in freeze-up and break-up dates, due to the above changes in atmospheric forcing, requires a temporal frequency of observations on the order of a few days, in order to capture the subsequent 1–2 week shifts in the timing of these events.

Studies such as this, which employ physically based models that include many of the relevant processes occurring in the natural system, can be used to guide the development of observational programs attempting to detect changes in complex, coupled land-atmosphere systems. In support of other studies, which have suggested that observations of freeze-up date, maximum ice thickness, and break-up date of lakes in middle- to high-latitude regions would be useful indicators of changes in local and regional climate, this study has begun to define the requirements of such an observational program.