1. Introduction
Lake ice growth and decay includes freeze-up in the autumn, a long period of growth and thickening in winter, a short period of ice melting and thinning and, finally, breakup and the complete disappearance of the ice cover in spring. The sensitivity of lake ice growth and decay, particularly freeze-up, break-up and thus duration (sometimes referred to as ice phenology), to climate variability and change has been demonstrated by observation on the ground and from space (e.g. Reference Palecki and BarryPalecki and Barry, 1986; Reference Robertson, Ragotzkie and MagnusonRobertson and others, 1992; Reference Wynne, Lillesand, Clayton and MagnusonWynne and others, 1998; Reference MagnusonMagnuson and others, 2000). The sensitivity of ice phenology and thickness to climate variability and change has also been investigated using computer models (e.g. Reference Liston and HallListon and Hall, 1995b; Reference Vavrus, Wynne and FoleyVavrus and others, 1996; Reference Stefan and FangStefan and Fang, 1997; Reference Zhang and JeffriesZhang and Jeffries, 2000). An important aspect of computer modelling of lake ice, indeed of any modelling exercise, is the evaluation of model performance by comparison of simulated results with recent observations and measurements (Reference Heron and WooHeron and Woo, 1994; Reference Liston and HallListon and Hall, 1995a; Reference Fang, Ellis and StefanFang and others, 1996; Reference Walsh, Vavrus, Foley, Fisher, Wynne and LenterWalsh and others, 1998; Reference Ménard, Duguay, Flato and RouseMénard and others, 2002; Reference Duguay, Flato, Jeffries, Ménard, Rouse and MorrisDuguay and others, 2003).
Using the Canadian Lake Ice Model (CLIMo; Reference Ménard, Duguay, Flato and RouseMénard and others, 2002; Duguay and others, 2003), Reference Morris, Jeffries and DuguayMorris and others (2005) simulate the effects of the mid-1970s Alaska climate shift on lake ice phenology, thickness and composition (congelation ice, snow ice) in central Alaska, USA, and the sensitivity of those variables to possible future changes in the magnitude and timing of air temperature and precipitation in the region. Prior to performing the climate sensitivity study, we assessed the performance of CLIMo by simulating ice phenology, thickness and composition in ice seasons 2001/02 and 2002/03 and comparing the results with observations and measurements made at ponds in the vicinity of Poker Flat Research Range (PFRR) and in Fairbanks, central Alaska. This paper describes the model performance assessment, which is based on a more comprehensive dataset than was available to Reference Duguay, Flato, Jeffries, Ménard, Rouse and MorrisDuguay and others (2003) for PFRR ice season 1999/2000. This assessment includes freeze-up and breakup data obtained every 1–2 days at 11 ponds, and weekly data for ice thickness and the depth and density of the snow on the ice during the winter.
2. Description of Model and Study Locations
CLIMo is an adaptation of the one-dimensional thermodynamic sea-ice model of Reference Flato and BrownFlato and Brown (1996). A comprehensive description of CLIMo is given in Reference Ménard, Duguay, Flato and RouseMénard and others (2002) and Reference Duguay, Flato, Jeffries, Ménard, Rouse and MorrisDuguay and others (2003). The following brief description is from those sources.
CLIMo is based on the one-dimensional unsteady heat-conduction equation, with penetrating solar radiation, of Reference Maykut and UntersteinerMaykut and Untersteiner (1971), i.e.
where ρ is density (kgm–3), C p is specific heat capacity (J kg–1 K–1), T is temperature (K), t is time (s), k is thermal conductivity (Wm–1 K–1), z is the vertical coordinate (positive downward), F sw is the downwelling shortwave radiative energy flux (Wm–2), I 0 is the fraction of the shortwave radiative flux that penetrates the surface, α is albedo and K is the bulk extinction coefficient for penetrating shortwave radiation.
The surface energy budget is then computed from
where F 0 is the net downward heat flux (Wm–2) absorbed at the surface, є is surface emissivity, σ is the Stefan–Boltzmann constant (5.67 × 10-8Wm–2 K–4), F lw is the downwelling longwave radiative energy flux (Wm–2), F lat is the downward latent heat flux (Wm–2) and F sens is the downward sensible heat flux (Wm–2).
Simulations are run with the following input variables: mean daily air temperature (˚C), snow depth (m), wind speed (m s–1), relative humidity (%) and cloud cover (tenths). Each of these five variables is available from the US National Weather Service (NWS) station at Fairbanks (NWS-FAI). Air-temperature and wind-speed data are also available at PFRR, 35km northeast of NWS-FAI.
The snow cover in the model is represented as a single slab without compaction and redistribution. Default snow density values can be used, or, when they are available, mean measured dry-snow and wet-snow density values can be specified for winter and spring, respectively. The other input variables are the time-step (1–24 d–1), the mixed layer depth (m) of the water below the ice, and site latitude. Model output values include mean daily total ice, snow-ice and congelation-ice thicknesses, freeze-up and break-up dates, duration and all components of the surface energy balance.
At PFRR, there are three main study sites: MST pond (65.13° N, 147.45° W), 31.6 mile pond (65.14° N, 147.44° W) and 33.5 mile pond (65.15° N, 147.38° W). MST pond is 0.7 km from 31.6 mile pond, and 2.5 km from 33.5 mile pond. The PFRR meteorological instruments are 1.5 km southwest of MST pond. In addition to the three main ponds, eight other nearby ponds located between 29.5 mile and 36.6 mile Steese Highway were visited on a daily basis in autumn and spring to monitor freeze-up and break-up, respectively. The Fairbanks study site is Aurora pond (64.85° N, 147.75° W), 13 km northeast of NWS-FAI. MST pond is a natural water body. The other PFRR ponds and Aurora pond are flooded gravel pits. The PFRR ponds have dimensions of the order of 100–200m and are up to 5m deep. Aurora pond is 200m across and up to 17m deep.
During the early stages of ice growth in autumn, total ice thickness was measured by drilling every 1–2 days. Once the ice was thick enough to walk on safely, heated wire ice-thickness gauges (Reference Ramseier and WeaverRamseier and Weaver, 1975) were installed and total ice thickness was measured weekly until the onset of melt. The depth and density of the snow on the ice were also measured during the weekly winter visits to each pond. Every 1–2 days during the melt season, the heated wire gauges were used to measure total ice thickness until it was no longer safe to be on the ice. Shortly before the onset of melt, the thickness of snow-ice layers was measured after dry holes had been drilled in the ice along 100 m long transects extending across the ponds. Snow ice forms at the top of the ice cover when the weight of snow depresses the ice surface below the water level, water flows up through cracks and floods the base of the snow, and the resultant slush freezes. Snow ice lies above congelation ice that grows as water freezes on the bottom of the ice cover.
3. Results
3.1. Ice thickness
PFRR simulations were run using a 2 m mixing depth, and cold- and warm-snow density values of 145 and 254 kgm–3, respectively (the values are the mean of all measurements at MST, 31.6 mile and 33.5 mile ponds). Aurora pond simulations were run using a 10m mixing depth, and cold- and warm-snow density values of 155 and 316 kgm–3, respectively.
For each location, the model was initially run with NWS-FAI data for the five meteorological input variables. The results (simulation 1) for the PFRR ponds and Aurora pond are shown in Figures 1 and 2, respectively. At each location, with the exception of the early stages of ice growth in October at PFRR and in mid-November at Aurora pond, there is poor agreement between simulated and measured ice-thickness values in the growth and decay seasons. At PFRR, the agreement between simulated and measured ice-thickness values in both the growth and decay seasons improves greatly when air temperature and wind speed measured at PFRR, and snow depth measured on the PFRR ponds, are substituted for NWS-FAI air-temperature, wind-speed and snow-depth input data (simulation 2, Fig. 1).
Simulated and measured ice thickness at MST pond, 31.6 mile pond and 33.5 mile pond, PFRR, October 2001–May 2002.
Simulated and measured ice thickness at Aurora pond, Fairbanks, November 2002–May 2003.
Simulation 2 agrees particularly well with the MST pond ice-thickness measurements (Fig. 1). Modifying the PFRR air-temperature input values by –2.5˚C and –4.5˚C results in simulations 3 and 4, which agree well with measured ice thickness at 33.5 mile and 31.6 mile ponds, respectively. We have no air-temperature data for the ponds, but the simulated differences are plausible, based on our experience of conditions during fieldwork visits. Unlike 31.6 and 33.5 mile ponds, which are surrounded on all sides by dense stands of mature, tall spruce, birch and aspen trees only 5–10m away from the water’s edge, MST pond is more open, with less dense stands of younger trees further away from the water’s edge. Consequently, MST pond receives more solar radiation than the other ponds, and is less sheltered from winds that mix the atmosphere and break down any temperature inversion immediately above the snow and ice.
In the case of Aurora pond, the primary factor behind the discrepancy between simulation 1 and the measurements (Fig. 2) is the representativeness of the NWS-FAI snow-depth data. Between the time of initial ice formation (8 November) and 12 December, the ice at Aurora pond was bare and uninsulated. Snow did not begin to accumulate on the ice until 13 December. When snow depth on the ice at Aurora pond is substituted for the NWS-FAI snow-depth data, the result is a significant improvement between simulated (simulation 2) and measured ice thickness. Using Aurora pond snow-depth data, the model simulates more rapid ice growth and thus greater ice thickness than when NWS-FAI snow-depth data are used.
3.2. Ice composition
Snow ice did not form at Aurora pond in 2002/03, and the model did not simulate any snow-ice formation. At the PFRR ponds, the thickest snow-ice layers occur at the margins of the ice cover, where snowdrifts accumulate and the weight of snow is thus greater than at the center of the ice cover. To avoid skewing the mean snow-ice layer thickness, we exclude the thick snow-ice layers due to edge effects and see that there was zero to negligible snow-ice formation across most of each pond in 2001/02 (Table 1). The model overestimates the amount of snow ice at MST pond, and correctly simulates the amount at 31.6 and 33.5 mile ponds (Table 1). Simulations of ice composition in other winters at PFRR have underestimated the amount by as much as 20%.
Mean measured and simulated thickness of snow ice at MST, 31.6 and 33.5 mile ponds, spring 2002. Superscripts 2–4 refer to simulations 2–4 in Figure 1
3.3. Freeze-up, break-up and duration
For PFRR simulations 2–4 (Fig. 1), the earliest and latest freeze-up dates are 14 and 18 October 2001, respectively. At the 11 ponds that were monitored, freeze-up occurred 5 days earlier, between 9 and 13 October. Freeze-up at the three main ponds covered this entire range. For both Aurora pond simulations (Fig. 2), freeze-up occurs on 13 November. Observed freeze-up was 5 days earlier, on 8 November.
For PFRR simulations 2–4, the earliest and latest break-up dates are 13 and 28 May 2002. At the 11 ponds, break-up occurred between 19 and 25 May. Break-up at the three main ponds occurred between 19 and 23 May. For Aurora pond simulations 1 and 2, break-up occurs on 30 April and 3 May, respectively. Observed break-up was on 10 May.
At PFRR, simulated break-up begins earlier and extends over a longer period of time than has been observed. At Aurora pond, simulated break-up is 7–10 days earlier than has been observed. Figure 3 suggests that the premature simulated break-up is related, in part, to the ice-thinning rate in spring. The model consistently overestimates the thinning rate by 10–20mmd–1. This is discussed further in section 4.
Simulated and measured daily thinning rates during the decay season April–May 2003 at (a) MST pond, 31.6 mile pond and 33.5 mile pond, PFRR, and (b) Aurora pond. The simulation in (a) was run using PFRR air-temperature, wind-speed and snow-depth data, and NWS-FAI cloud-cover and relative humidity data. Numbers in parentheses in the legend indicate the average amount by which the model overestimates the thinning rate at each location.
For PFRR simulations 2–4, ice duration varies between 211 and 224 days, while ice duration observed at the 11 ponds varied between 218 and 225 days. The model underestimates ice duration at MST, 31.6 and 33.5 mile ponds by 1.8%, 0.4% and 1.8%, respectively, relative to simulations 2, 3 and 4. For Aurora pond simulations 1 and 2, ice duration is 169 and 172 days, respectively. Observed duration was 183 days. Simulations 1 and 2 underestimate ice duration at Aurora pond by 8.3% and 6.4%, respectively.
4. Discussion
The model generally simulates a late freeze-up compared to observations. This is probably because the rate of heat loss from the water is too slow and ice formation is delayed. Heat loss from the water is determined by the open-water surface energy balance, which contains a user-designated mixed layer depth that is a measure of the heat content of the water body (Reference Duguay, Flato, Jeffries, Ménard, Rouse and MorrisDuguay and others, 2003). Fixing a mixed layer depth is a relatively simple approach compared to other lake sub-models which simulate not only the surface energy balance but also the water temperature profile and mixing (e.g. Reference Liston and HallListon and Hall, 1995a; Reference Fang, Ellis and StefanFang and others, 1996; Reference Vavrus, Wynne and FoleyVavrus and others, 1996).
Once the model has formed ice, the agreement between simulated and measured ice thickness is a function of the source of data used as model inputs. Using only NWS data to run the model yields poor agreement between simulations and measurements. This is primarily a matter of the representativeness of those data. For example, from early October to mid-April, PFRR is as much as 10˚C cooler than NWS-FAI, which is subject to an urban heat-island effect (Reference Magee, Curtis and WendlerMagee and others, 1999). Running the model with available PFRR data greatly improves the accuracy of the simulations. Located in a residential–industrial area 2.2 km from the city centre, Aurora pond is subject to the heat-island effect. Aurora pond ice-thickness simulations are sensitive only to the choice of snow input file.
During the melt season, the model overestimates the thinning rate when the best available locally relevant air-temperature data are used as inputs. This suggests that other factors play a role in the premature simulated breakup of the ice cover. We propose that the ice albedo is one important factor. CLIMo parameterizes albedo as a function of ice thickness (Reference Duguay, Flato, Jeffries, Ménard, Rouse and MorrisDuguay and others, 2003), thus ignoring the composition of the ice. If we calculate albedo values from the simulated shortwave radiation reflected off the melting ice surface and the total incoming shortwave radiation as determined by the cloud-cover input data, we obtain mean simulated albedo values of 0.17±0.02 (min. 0.15, max. 0.22) and 0.16± 0.01 (min. 0.15, max. 0.18) for the periods illustrated in Figure 3a and b, respectively.
The simulated albedo values are lower than those measured for snow ice and congelation ice. Snow ice has an albedo of ~0.4 (Reference BolsengaBolsenga, 1983) and retards the thinning rate during the early decay season (Reference AshtonAshton, 1986, p. 355). The albedo of congelation ice varies according to its crystal texture. Melting congelation ice with a columnar texture has an albedo of 0.2–0.4 (Reference Prowse and MarshProwse and Marsh, 1989; Reference Heron and WooHeron and Woo, 1994), and that with a non-columnar texture an albedo of 0.4–0.55 (Heron and Woo, 1994). Both types of congelation ice, often overlain by snow ice, occur at PFRR and Aurora pond. As a consequence of its higher albedo, the ice thins more slowly and break-up is delayed relative to the simulations.
A better knowledge of melting lake ice albedo and improved albedo parameterizations would likely improve the simulation of ice-thinning rate and break-up in spring. Here we have seen that simulated break-up occurs too early compared to observations because of accelerated thinning. Coupled with freeze-up being simulated late compared to observations, the simulated ice duration is too short. However, if ice duration is examined in relative terms rather than absolute number of days, the model actually underestimates ice duration only by <10%.
This simulation of ice duration is actually quite good in view of the fact that the definitions of simulated and observed freeze-up and break-up are not the same. Simulated freeze-up and break-up occur in a one-dimensional space, and are defined as the first day of ice formation and first day of zero ice thickness, respectively. Observations of freeze-up and break-up are made in three-dimensional space, and refer to the time when the water body has a 100% ice cover and zero ice cover, respectively. As we have observed at PFRR and Aurora pond, it often requires a few days for even a small water body to achieve 100% ice cover in autumn as the area of ice becomes progressively larger. Likewise, in spring the area of ice decreases as the area of open water increases until it reaches 100%. Simulated freeze-up and break-up are instantaneous events defined by ice thickness. Observed freeze-up and break-up are progressive events defined by area. Consequently, a model is unlikely to achieve perfect agreement with observations of freeze-up, break-up and duration.
For winter 2001/02, the model achieved good results with respect to ice formation and composition. Simulations for other winters have tended to underestimate the amount of snow ice. Comparisons of the simulated and observed ice composition are also affected by the dimensions of the model and the observations. Snow ice often does not form a uniform layer across a lake ice sheet (e.g. Reference Adams and ProwseAdams and Prowse, 1981), and a one-dimensional model cannot reproduce the resultant high spatial variability of snow-ice thickness and its proportion of the total ice thickness. Multiple measurements of snow-ice layer thickness across a lake ice sheet, such as we make at Poker Flat, can be reduced to a single value (e.g. the average thickness), but there is a high probability that the single snow-ice quantity provided by a one-dimensional model will not agree with this amount. The lack of agreement between simulated and observed snow ice can also be attributed to model factors. The absence of compaction during snow accumulation and capillary action during flooding (i.e. the model underestimates the depth of wetted snow that is transformed into snow ice) also explains the underestimation of snow ice by CLIMo (Duguay and others, 2003).
Conclusion
CLIMo simulates ice growth and decay in a complex air–snow-ice–water system. This study, and those of Reference Ménard, Duguay, Flato and RouseMénard and others (2002) and Reference Duguay, Flato, Jeffries, Ménard, Rouse and MorrisDuguay and others (2003), show that overall the model does this very well. The best agreement between simulations and observations is achieved when the model is run with input data that represent as closely as possible the environmental conditions at the ponds being investigated. The level of agreement between the simulations and field data is sufficient to conclude that, where good meteorological and other data exist to run the model, it can be relied on to make accurate predictions of past variability of lake ice growth and decay, and the response of freeze-up, ice thickness and composition, break-up and duration to future climate change.
There are discrepancies between simulations and observations, which suggest that improvements could be made to the lake, snow and snow-ice sub-models, and the melting-ice albedo parameterization. Data on melting lake ice albedo are actually more sparse than the values presented in section 4 might suggest; systematic, time-dependent measurements, and experimentation with albedo parameterizations would be a valuable contribution to lake ice studies. Some discrepancies between simulations and observations are unavoidable when one-dimensional model results are being compared with three-dimensional observations. Further progress in lake ice geophysics would benefit from two- and three-dimensional modelling of ice processes. This would allow factors such as lake depth and area, which can have a significant effect on freeze-up and break up, to be taken into account and improve simulations. Coupled with physical and ecological models of lakes, such ice models would add to the knowledge and understanding of variability and change in frozen lake interactions and feedbacks.
Acknowledgements
The Poker Flat results presented in this paper were made possible with the support of US National Science Foundation (NSF) grant OPP 0117645 and the International Arctic Research Center (IARC) at the University of Alaska Fairbanks (NSF Cooperative Agreement OPP 0002239). We thank the Bennett and Irving families for making the measurements at Aurora pond under the auspices of the Alaska Lake Ice and Snow Observatory Network, a program supported by the University of Alaska Natural Resources Fund, NSF grant OPP 0326631 and IARC (NSF Cooperative Agreement OPP-0002239). CLIMo was developed with support from the Meteorological Service of Canada. We thank the three anonymous reviewers and the Scientific Editor, D. MacAyeal, for constructive reviews that helped us complete a concise and focused paper.
