An energy-balance model of lake-ice evolution

Abstract

A physically based mathematical model of the coupled lake, lake ice, snow and atmosphere system is developed for studying terrestrial-atmospheric interactions in high-elevation and high-latitude regions. The ability to model lake-ice freeze-up, break-up, total ice thickness and ice type offers the potential to describe the effects of climate change in these regions. Model output is validated against lake-ice observations made during the winter of 1992–93 in Glacier National Park, Montana. U.S.A. The model is driven with observed daily atmospheric forcing of precipitation, wind speed and air temperature. In addition to simulating complete energy-balance components over the annual cycle, model output includes ice freeze-up and break-up dates, and the end-of-season clear ice, snow-ice and total ice depths for two nearby lakes in Glacier National Park, each in a different topographic setting. Modeled ice features are found to agree closely with the lake-ice observations.

Model simulations illustrate the key role that the wind component of the local climatic regime plays on the growth and decay of lake ice. The wind speed affects both the surface temperature and the accumulation of snow on the lake-ice surface. Higher snow accumulations on the lake ice depress the ice surface below the water line, causing the snow to become saturated and leading to the formation of snow-ice deposits. In regions having higher wind speeds, significantly less snow accumulates on the lake-ice surface, thus limiting snow-ice formation. The ice produced by these two different mechanisms has distinctly different optical and radiative properties, and affects the monitoring of frozen lakes using remote-sensing techniques.

Introduction

Evidence from the available air-temperature record indicates a distinct two-decade warming trend over northern land areas during the winter and Spring seasons (Chapman and Walsh, 1993). In addition to conventional meteorological observations, monitoring and analysis of lake-ice parameters such as freeze-up date, maximum ice thickness and break-up date have been proposed as useful indicators of regional climate change (Palecki and Barry, 1986; Robertson and Others, 1992). Since frozen lakes provide an index of integrated seasonal temperature trends, the analysis 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.

While numerous previous researchers have studied the relationship between air temperature and lake-ice formation and break-up, a quantitative understanding of the underlying physical processes which lead to the timing of these climatic events is limited. Researchers have adopted various means to quantify the relationships between air temperature and lake-ice cover. Rannie (1983) and palecki and Barry (1986) performed regression analyses between key atmospheric and ice-related Variables such as air temperature, freeze-up dale, break-up date and ice duration for lakes where significant meteorological and ice-cover data sets exist. Other studies have used an integrated average air-temperature index to predict freeze-up and break-up dates (McFadden, 1965). In such an approach, the dales of freeze-up and break-up are assumed to occur when the index reaches a specified threshold.

A more complex approach to describing lake-ice processes is through the formulation of physically based models which account for the relevant process occurring within, and at the boundaries of the lake, ice, snow and atmospheric components of the natural system. Following the original sea-ice research Conducted by Stefan in the late 1800s, mathematical descriptions of freezing and thawing processes are often referred to as Stefan problems. In previous lake-ice studies, Rodhe (1952) and Bilello (1964) modeled lake-ice formation by considering the sensible-heat exchange at the water air interface. Robertson and others (1992) adopted a similar approach to relate lake freeze-up and break-up dates to air temperature.

Maykut and Untersteiner (1971) developed a thermodynamic model which included an accounting for the energy fluxes within snow and ice layers, and at the snow air interface. Their model was used to describe the seasonal evolution of sea ice. Patterson and Hamblin (1988) presented a two-component (snow and ice) thermodynamic model based on the formulation of Maykut and Untersteiner (1971). They coupled the snow- and ice-component sub-model to a reservoir-mixing model, thus providing an interactive linkage between the ice cover and the underlying body of water. Gu and Stefan (1990) used a snow and ice model similar to that of Patterson and Hamblin (1988). Their lake model simulates the year-round hydrothermal processes and includes sediment heat fluxes occurring at the lake bottom. Marsh (1991) used a simple lake-ice-growth model based on a modified Stefan’s equation. Various solutions to the Stefan problem of melting and freezing have been discussed by Lock (1990). Using an energy-balance model, Heron and Woo (1994) studied melt processes during the decay of a lake-ice cover in the Arctic. General discussions of the thermal regime of ice-covered lakes can be found in Pivovarov (1973) and Ashton (1986).

In this paper, an energy-balance model describing the seasonal evolution of lake ice is developed and is used to further the understanding of the interactive role that the atmosphere and the unfrozen water body plays in determining ice-growth mechanisms and other energy-related processes. Model output is validated against lake-ice observations made during the winter of 1992–93 in Glacier National Park, northern Montana, U.S.A. In addition, model output is used to assist in explaining variations in synthetic aperture radar (SAR) signatures obtained from several lakes in Glacier National Park during the same period (Hall and others, 1994). The sensitivity of this lake-ice model to variations in atmospheric forcing has been described by Liston and Hall (1995).

The one-dimensional, unsteady model is composed of four major sub-models. First, a surface-energy balance sub-model is implemented to determine the surface temperature and energy available for freezing or melting. Secondly, 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. Thirdly, a snow sub-model describes the snow depth and density as it accumulates, metamorphoses and melts on top of the lake ice. Fourthly, a lake-ice growth sub-model produces ice by two different mechanisms: (a) relatively bubble-free clear ice grows at the ice-water interface due to thermal gradients within the ice, and (b) 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. In addition to complete energy-balance components over the annual cycle, key model output includes the dates of ice freeze-up and break-up, and the end-of-season clear-ice, snow-ice and total ice depths. The model is driven with observed daily atmospheric forcing of precipitation, wind speed and air temperature. A schematic illustrating the components of the model and a representative temperature profile is given in Figure 1).

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

2. Surface-Energy Balance

The surface energy-balance equation is

(1)

where Q_si is the solar radiation reaching the surface of the Earth, Q_li is the incoming longwave radiation, Q_le is the emitted longwave radiation, Q_h is the turbulent exchange of sensible heat, Q_c 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. All of the energy terms, Q, have units of W m⁻².

When coupled to the coupled lake, lake-ice and snow models, the solution of this energy balance prov ides 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. In addition to the following equations describing the relevant physical processes, observed daily atmospheric forcing of wind speed, precipitation and maximum and minimum surface-air temperature are required to solve the energy balance.

The solar radiation striking a horizontal surface, Q _si, is given by

(2)

where S* is the solar irradiance at the top of the atmosphere striking a surface normal to the solar beam (= 1370 W m⁻²) (Kyle and others, 1985) and γ is the net sky transmissivity or the fraction of solar radiation that arrives at the surface. The solar elevation angle, α, is defined as

(3)

where ϕ is latitude, T is the hour angle measured from local solar noon and δ is the solar-declination angle approximated by

(4)

where ϕ is the latitude of the Tropic of Cancer (23.45°N), d is the Gregorian Day, d _r is the day of the summer solstice (173) and d _y is the average number of days in a year (365.25).

To account for the scattering, absorption and reflection of shortwave radiation by clouds, the solar radiation is scaled according to

(5)

where σ_c, represents the fraction of cloud cover (Burridge and Gadd, 1974). Because of lack of available observations, in the current application, a constant cloud-coverage fraction of 0.4 has been assumed (Conway and Liston, 1974). This assumption of constant cloud fraction was relaxed in Liston and Hall (1995) during simulations of ice growth and decay on Back Bay, Great Slave Lake, northern Ganada.

An analytical downward longwave-radiation equation which considers clear skies and standard atmospheric conditions has been developed by Brutsaert (1975). To correct apparent deficiencies in this formulation at air temperatures below 0°C (Aase and Idso, 1978), Satterlund (1979) suggested the empirical description

(6)

where T _a, is the surface-air temperature and e _a is the atmospheric vapor pressure. In the current study, no attempt has been made to modify this formulation for Q _li to account for the presence of clouds. The average daily surface-air temperature, T _a, is simply assumed to be the arithmetic mean of the daily maximum and minimum temperatures, T_max and T_min, respectively.

In an analysis of observed dew-point and night-minimum temperatures, Hungerford and others (1989) concluded that in Montana the daily minimum temperature is a reasonable approximation of the dew-point temperature. (In locations where dew-point temperatures are not frequently reached on a daily basis, due, for example to arid conditions or minimal diurnal temperature fluctuations, this approximation would not be valid. Adopting this assumption, the atmospheric vapor pressure. e_a, is computed from the minimum temperature using the formula (Fleagle and Businger, 1980)

(7)

The longwave radiation emitted by the snow surface is computed under the assumption that snow emits as a grey body

(8)

where T _s0 is the snow-surface temperature, σ is the Stefan-Boltzmann constant and σ_s is the snow emissivity, assumed to be 0.98.

The turbulent exchange of sensible and latent heat, Q _li, and Q _e, respectively, have been given by (Price and Dunne, 1976)

(9)

(10)

where ρ_a is the density of the air, T _a is the surface-air temperature, e _a is the vapor pressure of the air, e _s0 is the vapor pressure of the surface, Ç is a non-dimensional stability function, ρ_a is atmospheric pressure, C _p is the specific heat of air and L _s is the latent heat of sublimation, D _h and D _e are exchange coeficients for sensible and latent heat, respectively

(11)

where k is von Karman’s constant, U _z is the wind speed at reference height z and z ₀ is the roughness length. The roughness length varies with surface type and in this application z ₀ is assumed to be 0.001, 0.005 and 0.000001 m over open water, snow and bare ice, respectively. under stable atmospheric conditions, the stability function ζ modifies the turbulent fluxes through the formula

(12)

where the Richardson number for the atmosphere, Ri_a is given by

(13)

g is the gravitational acceleration and the atmospheric temperature gradient is computed using the reference-level air temperature and the surface temperature, T _s0.

A representative atmospheric pressure, P_a, for the location of interest is given by

(14)

where P ₀ is a reference sea-level pressure (101 300 Pa), H is the scale height of atmosphere (≈8000 m) and Z _g is the elevation of the study site (Wallace and Hobbs, 1977).

Heat conducted through the ice and snowpack is

(15)

where k is the effective thermal conductivity of the snow-ice-water matrix. In the model, for the case where lake ice is present, the conductive heat flux is approximated by

(16)

where T _f is the water-freezing temperature, z is the depth of each layer and the subscripts i, s, m and w indicate individual layers of ice, snow, a mix of snow and water, and water, respectively. For the case of a lake with no ice, the flux contribution due to conduction is

(17)

where the subscript 1k indicates a lake value. In the model, conditions within the top 1 m of the lake are used to compute the conduction. The determination of many of these conduction parameters requires a snow model, a lake-ice modal and an ice-free lake model. These sub-models are described later in this paper.

The albedo, α_s, decreases with increasing snow density, ρ_s, following Anderson (1976)

(18)

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

(19)

Several researchers have indicated that there appears to be no dependency of snow albedo on snow density (Winther, 1993). While in addition, the decrease in snow albedo with time since the last snowfall, corresponding to an increase in snow grain-size, is also well noted. In our formulation, increases in snow density with time are assumed to correspond to increases in grain-size, which lead to changes in albedo according to the preceding formulation. The evolution of the snow density is described in the following section. Lake ice is assumed to have an albedo of 0.25 and a lake surface without snow or ice has an albedo of 0.06.

To implement the surface energy-balance model, daily values of T _max, T _min and U _z are assumed to be provided from observations or some suitable atmospheric model. These data, in combination with the preceding energy-balance equations and the snow, lake-ice and lake models described below, allow computation of the surface temperature. To solve the system of equations for the surface temperature, the equations are east in the form

(20)

and solved iteratively for T _s0 using the Newton-Raphson method.

In the presence of snow and/or ice, surface temperatures T < 0°C resulting from the surface-energy balance indicates 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 _f, liquid water present on the lake-ice surface or within the snowpack,.

3. Lake-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 ice-water 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

(21)

Where ρ_i is the ice density, h _w is the convective-transfer coefficient, assumed to be 0.56 W K⁻¹ m⁻² for a calm lake based on conductivity considerations, T _f is the water-freezing temperature, T _W is the lake-water temperature at a depth of 0.33 m as determined from the lake model described below, Z _i is the ice depth, Z _s is the snow depth, K _i is the thermal conductivity of ice, and

is the velocity of the moving ice boundary. For the largely snow-cowed lakes considered in this study, the influence of solar radiation penetrating the ice has been assumed to be of secondary importance in this formulation.

Solution of the previous 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,

(22)

where this relationship has been fitted to a collection of data by several researchers presented by Mellor (1977) (Verseghy, 1991).

Density changes occur in the snowpack by two mechanisms. First, density increase can result from compaction and, following Anderson (1976),

(23)

where T _s is the snow temperature, W _s is the weight of snow above a given layer expressed in water-equivalent depth, and Α ₁ and A ₂ are constants set equal to 0.0013 m⁻¹ s⁻¹ and 0.021 m³ kg⁻¹, respectively, based On Kojima (1967). In the current application, the snow is considered as a single layer and the snow temperature is determined from

(24)

where the temperature at the snow-ice interface, T _si, is given by

(25)

as determined from energy-balance considerations under steady-state conditions. For this single-layer implementation, the weight of die snow, W _s, is defined as half of the snow water equivalent, W, where

(26)

and ρ_w is the density of water.

The second density-modifying process results from snow melting. The melted snow reduces the snow depth and is redistributed through the snowpack until a maximum snow density, assumed to be 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 potential snow melt, M _p, given by

(27)

is used first to melt available snow and then the lake ice. The energy available to freeze water within the snowpack, or on the ice surface, is used to compute the potential freeze depth, F _p,

(28)

which is used to add to the lake-ice depth either through freezing water held within the slush mixture of snow and water, or through the freezing of water existing on the lake surface in the absence of snow. Throughout the melting and freezing processes, snow, ice and water depths, and snow density, are adjusted accordingly, while taking into account the density differences between the solid and liquid phases.

Given input of liquid-equivalent precipitation, the precipitation is assumed to fall as snow if the wet-bulb temperature, T _wb, is lower than 1°C, where the wet-bulb temperature is given by (Rogers, 1979)

(29)

This equation is solved iteratively for T _wb, using the Newton Raphson method.

Precipitation falling as rain contributes to the liquid-water store on the lake surface. Snow falling on a frozen lake surface supporting more water than that required to saturate the existing snow cover adds to the layer of surface water. Snow falling an a bare-ice surface or an unsaturated snow surface accumulates as new snow. Following Anderson (1976), based on data by LaChapelle (1969), the new snow density, ρ_us, is given by

(30)

This new snow is added to the existing snowpack where 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 completely the overburden of snow cover at a level where the top of the ice equals that of the water. Buoyancy considerations indicate that 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

(31)

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 are assumed to form or be present, which allow lake water to saturate the ice surface. Under this condition, the solid fraction of the snow cover adds to the upward-buoyancy force and the water which had previously accumulated on the ice surface is inter-connected with the lake water. This leads to a new force balance

(32)

where z _m, is the depth of the water-saturated slush mixture. This can be solved for the depth of the slush-layer setting above the lake-ice surface

(33)

The liquid-water store, z _w, within this slush layer is given by

(34)

and is available to form snow-ice depending on the available freezing energy, Qf.

4. 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

(35)

where Τ _w is the water temperature, t is time, z is the vertical coordinate measured downward from the water surface, K _H is the turbulent diffusivity 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

(36)

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

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

(37)

where the Richardson number for the lake water, Ri_w, is

(38)

the surface-shear velocity, w _s ^*, is

(39)

and the shear-velocity decay constant, K*, is given by

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

The Brunt-Vaisala frequency, N ⁻², is determined by the density stratification

(41)

Computation of N ⁻² requires a description of the water-density variation with temperature. Kell (1972) presented a best fit of measurements for water for the temperature range 0–150°C in the form of a seven-parameter function

(42)

where the a and b constants can be found in the Kell (1972) reference. 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. The influence of model spin-up is accounted for by starting the model early enough during the summer/fall that the lake comes into balance with the atmospheric forcing and the freeze-up date is unaffected by the model-initiation date. A couple of months lead time was found to be more than adequate. In the absence of ice, the lake surface-temperature boundary condition is determined from the surface-energy balance; if ice is present, the temperature of the ice-water interface equals the water-freezing temperature. The heat flux from lake-bottom sediments is recognized to play an important role in governing the thermal regime of small, shallow lakes. For the relatively large and deep lakes addressed in this paper, a zero-flux boundary condition is appropriate and employed for the lower boundary

(43)

where the subscript z_max is the depth of the deepest lake-model layer, In the model, this boundary condition is applied at a depth of 42 m. The 42 m deep water column is divided into 25 individual layers ranging in depth from 0.167 m at the surface to 3 m at the bottom. The heat-transport equation is solved using a fully implicit, control volume-based, finite-difference scheme.

To initiate ice on the lake surface, a 1 mm thick ice layer is formed when the water temperature at a depth of 0.33 m is less than the water-freezing temperature. In the presence of lake ice, the surface-shear stress at the ice-water interface is zero and heat transfer within the lake below the ice cover occurs through molecular diffusive processes only; til this case K _H is zero.

5. Model Simulations

Two lakes in eastern Glacier National Park (GNP), Montana, are studied using the coupled lake-atmosphere model described in the preceding sections. Field observations of lake-ice depth and ice type (clear ice and snow-ice) were obtained for St. Mary and Lower Two Medicine Lakes during the winter of 1992–93. Lower Two Medicine Lake is located approximately 30 km south-southeast of St. Mary Lake (Fig. 2). Daily meteorological observations were collected at the St. Mary GNP ranger station which is located near the mouth of St. Mary Lake (Fig. 2). These meteorological data include daily maximum and minimum temperature and water-equivalent precipitation. In addition, observations of wind speed and direction were made daily at 1630 h local time. The meteorological data cover the period 1 September 1992 through 31 May 1993. The period of model integration coincides with these available atmospheric forcing data.

Fig. 2. The relationship between the St. Mary and Lower Two Medicine Lakes and the surrounding topography. The prevailing storm winds in this region arrive from the southwest.

The relationships between the two lakes and the surrounding topography are illustrated in Figure 2. St. Mary Lake is at an elevation of 1367 m with approximate dimensions of 1 km by 15 km. The upper end of the lake is aligned from West-southwest to east-northeast and lies at the foot of mountains which rise to as high as 1500 m above its surface to the north and south. To the west of the upper end of the lake lies a broad valley which leads to the Continental Divide. The lower end of the lake runs from southwest to northeast. The upper end of this lower part of the lake roughly aligns with a southwest-oriented valley which extends towards a relatively low pass through the Continental Divide. To the east of the lake the topography comprises gently rolling hills.

Lower Two Medicine Lake is at an elevation of 1488 m with approximate dimensions of 0.75 km by 5 km. The long axis of the lake is aligned northwest to southeast and lies at the foot of mountains which rise to as high as 1500 m above the lake surface on all sides except to the east where the topography is composed of gently rolling hills.

In this region of the northwestern United States the prevailing storm winds arrive from the southwest (Bryson and Hare, 1974; Wallace and Hobbs, 1977) (Fig. 2). As a consequence, the storm winds coming over the Continental Divide are directed down the southwest-aligned valleys which feed into St. Mary Lake, resulting in relatively strong winds on St. Mary Lake during storm events. In contrast. Lower Two Medicine Lake is surrounded by mountains to the south and west (Fig. 2). These mountains protect the lake from the prevailing storm winds and lead to much lower wind speeds over Lower Two Medicine Lake than over St. Mary Lake.

In the model integrations, the temperature data obtained from the St. Mary ranger station are assumed to be valid for both lakes. In addition, the observed wind speeds are applied directly for the St. Mary Lake situation. To account for the effect of surrounding topography on winds over Lower Two Medicine Lake, the wind speed measured at St. Mary is used after being multiplied by a factor of 0.4; wind-speed observations for Lower Two Medicine Lake are not available. This value has been chosen to provide the best agreement between the model predictions and the ice observations while still being consistent with the influence of the surrounding topography.

Snow accumulation on a relatively flat surface, such as that of a lake, is dependent upon complex inter-relationships between wind speed, air temperature and precipitation, in combination with historical factors which serve to define the shear strength of the existing snowpack. Unfortunately, models describing these relationships do not exist. To account for the obvious connection between reduced snow accumulation, over a flat lake, in response to higher wind speeds, the precipitation accumulating on the frozen lake surface is determined by multiplying the observed snow water-equivalent precipitation by the empirical scaling factor β, where β has been chosen to equal 0.1 and 0.9 for St. Mary and Lower Two Medicine Lakes, respectively. These values were chosen to optimize the agreement between model output and observations. An alternative to introducing this scaling factor is to prescribe the snow depth and density on the lake based on observations; this approach requires more frequent snow observations than are available for the current study. In summary, the only difference between the atmospheric forcing for the St. Mary Lake and Lower Two Medicine Lake integrations is that Lower Two Medicine Lake is assumed to have lower wind speeds and higher snow accumulations than those occurring on St. Mary Lake.

The St. Mary observations of temperature, wind speed and precipitation are found in Figure 3a- c, respectively. During October and November, the average daily air temperature reaches below-freezing temperatures on several occasions and stays below freezing through December and during almost all of January. In late January and early February, air temperatures rise above freezing before dropping during the last half of February. March begins with above-freezing temperatures and, after a cooler period in the middle of the month, the air temperature typically stays above freezing through May (Fig. 3a). The wind speed is highly variable throughout the period and ranges from near zero to 22 m s⁻¹ (Fig. 3b.) Precipitation events occur frequently and reach a maximum of 16 mm d⁻¹ (Fig. 3c).

Fig. 3. Measurements at St. Maty, Glacier National Park, Montana, U.S.A. a. Average daily temperature computed from the daily maximum and minimum air temperatures. b. Daily wind speed. Observations made at 1630 h local time. c. Daily snow water-equivalent precipitation.

The total ice depth simulated by the model for St. Mary Lake is compared with observations in Figure 4. The average of the observed depth values is included along with the maximum and minimum observations for each observation time. The number of observations varies from one field 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 of November. Ice break-up is modeled to occur during the first week of April. The high variability among depth observations for a given time can be partially explained by the variable depths of snowdrifts accumulating on the ice surface. 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 60 cm in late February. The relatively warm periods in early February and March, indicated by the temperature record, have led to the reductions in ice depth shown in Figure 4. The final melting of lake ice begins in the middle of March and continues until the ice is entirely melted in early April.

Fig. 4. Comparison between the total ice depth simulated by the model and the observations far St. Mary Lake. The average and range of the observed values are also indicated.

On St. Mary Lake, no snow-ice is produced by the model (Fig. 5). This curve is coincident with the zero line in the plot. Although the observations of snow-ice depth are quite variable, the madel-produced result lies within the observed values, also shown in Figure 5. In this figure the total ice depth from Figure 4 has been included for reference. Under conditions of minimal snow accumul-ation such as that represented by this St. Mary Lake simulation, the snow cover is rarely sufficient to depress the ice surface below the water line; if this should occur, it can happen only early in the freeze-up period when the ice is relatively thin. Consequently, snow-ice formation is expected to be limited in relatively cold and windy environments. The clear ice produced by the model is given by the difference between the total ice depth and the snow-ice depth. In the St. Mary Lake simulation, all of the ice produced by the model is clear ice.

Fig. 5. Comparison between the snow-ice depth simulated by the model and the observations for St.. Mary Lake. The average and range of the observed values are also indicated, along with a reference curve showing the total modeled ice depth obtained from Figure 4. (Note that the model did not simulate any snow-ice on St. Mary Lake.)

The model simulation for Lower Two Medicine Lake also agrees well with the measured total ice depth (Fig. 6). The number of measurements varies from one field excursion to the next and ranges from one to three observations per trip to the lake. Measurements were not always made in the same part of the lake during each trip. The modeled ice freeze-up date lags the observations by approximately 1 week. While the ice growth during December is slightly less than the model simulation, the general pattern of the measured ice growth is captured by the model. The remainder of the simulation continues to. follow the observations. In response to the different atmospheric forcing, which allowed the introduction of snow-ice growth not present on St. Mary Lake, the model produces a maximum total ice depth of 70 cm, 10 cm greater than that simulated for St. Mary Lake. The final continuous melting event begins in the middle of March and proceeds until the middle of April when all of the ice has gone. The lack of April ice observations precludes model validation during this period of the simulation.

Fig. 6. Comparison between the total ice depth simulated by the model and the observations for Lower Two Medicine Lake. The average and range of the observed values are also indicated.

In contrast to St. Mary Lake, the ice measured on Lower Two Medicine Lake is composed of roughly 50% snow-ice (Fig. 7). Also included in Figure 7 are the snow-ice observations and the total ice depth obtained from Figure 6. The clear-ice depth is given by the difference between the total ice and snow-ice depths. Again, the modeled snow-ice profile is within the variability of the observations. The increased snow accumulation in response to the lower wind speeds at Lower Two Medicine Lake has initiated the snow-ice growth mechanism. Snow accumulation an the lake surface has been sufficient to depress the ice surface below the lake-water line, thus saturating some fraction of the overlying snow cover. This slush has then frozen, thus producing the snow-ice layers. The reductions in ice depth due to melting in early February and March, which are evident in the St. Mary Lake simulation (Fig. 4), are not found in the Lower Two Medicine Lake simulation (Fig. 6). In this case, the melting energy has gone towards first melting the snow cover present on the Lower Two Medicine Lake surface, and the surplus energy was insufficient to reduce the ice depth appreciably.

Fig. 7. Comparison between the snow-ice depth simulated by the model and the observations for Lower Two Medicine Lake. The average and range of the observed values are also indicated, along with a reference curve showing the total modeled ice depth obtained from Figure 6.

The differences between the ice-growth mechanisms in the two lakes produce ice of distinctly different optical and radiative properties. Because of the abundance of small bubbles in the snow-ice, this ice form is opaque and frequently referred to as white ice. In contrast, the ice produced by thermal growth is clear and contains relatively few bubbles. Studies have shown that the density of snow-ice is often within a few per cent of the clear-ice density (Ager, 1962).

Lake-ice studies in northern Alaska using synthetic aperture radar (SAR) have suggested that the within-ice bubble structure is partially responsible for variations in radar back-scatter (Weeks and others, 1978, 1981; Wakabayashi and others, 1993). As part of a companion study to the modeling work described in this paper, SAR data covering Glacier National Park have been obtained for the months of November 1992 through March 1993 (Hall and others, 1994). The distinctly different bubble structures found within the lakes of Glacier National Park are thought to explain, at least in part, the variations between the radar signatures obtained from various lakes within the park.

These model simulations of Montana lakes have shown that the formation and growth of lake ice is sensitive to wind speed, which is strongly dependent upon local topography. Under both light and strong winds, lakes were found to produce ice through thermal growth at the ice-water interface but, under conditions of light winds, the formation of snow-ice is enhanced. Higher snow accumulations on the lake-ice surface produce a condition where the ice surface is depressed below the water line, causing the snow to become saturated and leading to the formation of relatively thick snow-ice deposits. In regions having wind speeds high enough to limit snow accumulation, and/or where precipitation is relatively low, the reduced snow accumulation on the lake-ice surface inhibits snow-ice formation.

6. Concluding Remarks

Physically based mathematical models of the coupled lake, lake-ice, snow and atmosphere system provide a mechanism by which to study terrestrial—atmosphere interactions for the climate-sensitive northern regions. With such a model, freeze-up and ice break-up dates can be predicted with reasonable accuracy, provided that the appropriate atmospheric forcing is known. Given observations of the freeze-up and ice break-up dates of specific lakes, information regarding the local and regional climatic regime of that lake can be inferred.

This study illustrates the key role that the wind component of the local climatic regime plays on the growth and decay of lake ice. The wind speed affects both the surface temperature, obtained from energy-balance considerations, and the accumulation of snow on the lake-ice surface. Snow-accumulation amount plays a significant rule in determining the relative importance of lake-ice growth through the thermal and snow-ice growth mechanisms. The ice produced by these two different mechanisms has distinctly different optical and radiative properties and is, consequently, thought to affect radar signatures (Hall and others, 1994).

Further research with physically based models are required to evaluate the sensitivity of lake freeze-up and ice break-up dates to long-term, sustained changes in winter climatic conditions. Such a study would assist in establishing the feasibility of tracking past and future climate trends by using lake Freeze-up and beak-up dates as an integrated index of fall and winter temperatures, and other- key climatic variables.