Wed, 09 Apr 2014 in Journal of Limnology
One-Dimensional Simulation of Lake and Ice Dynamics During Winter
Abstract
An ice formation model, based on the solution of the heat conduction equation across blue ice, white ice and snow cover, is integrated into the Dynamic Reservoir Simulation Model (DYRESM) to allow for one-dimensional (vertical) winter simulation of lake dynamics during periods of ice cover. This is an extension of a previous three-layer snow and ice model to include two-way coupling between the ice and the water column. The process-based ice formation is suitable for application to mid-latitude regions and includes: snowmelt due to rain; formation of white ice; and variability in snow density, snow conductivity, and ice and snow albedo. The model was validated against published observations from Harmon Lake, British Columbia, and new observations from Eagle Lake, Ontario. The ice thickness and water column temperature profile beneath the ice were predicted with Root Mean Square Deviations (RMSD) of 1 cm and 0.38°C, respectively, during the winter of1990-91in Harmon Lake. In Eagle Lake the 2011-12 year-round water column temperature profile was predicted with an RMSD of 1.8°C. Improved prediction of under-ice lake temperature, relative to published results from simpler models, demonstrates the need for models that accurately capture ice-formation processes, including ice to water column coupling, formation of both blue and white ice layers, and process-based ice and snow parameters (density, conductivity and albedo).
Main Text
INTRODUCTION
To numerically model lake dynamics in winter, there is a need for the development of ice-cover algorithms. These algorithms must accurately reproduce the capacity of ice cover to reduce the surface wind stress and both thermodynamic and biogeochemical surface fluxes. Consequently, the hydrodynamic processes under ice have slower characteristic time-scales than in open water; however, ice-covered lakes are not completely quiescent (Kirillin et al., 2012). Neglecting these ice-effects on free-surface momentum and thermodynamic fluxes results in over-prediction of the spring surface temperature (Shore, 2009; Wang et al., 2010) and excessive mixing during the period of ice cover (Farmer, 1975; Kirillin et al., 2012).
Three-dimensional (3D) Reynolds-averaged lake models have recently implemented ice-cover algorithms (Wang et al., 2010; Oveisy et al., 2012), that are suited for seasonal process-oriented simulations. However, for long-term studies (tens to hundreds of years) involving coupling between climate and lake biogeochemistry, 3D models are computationally prohibitive and one-dimensional (1D) vertical models with ice cover are required (see discussion in Sahoo and Schladow, 2008; Mackay et al., 2009; Yang et al., 2012). Several such 1D models have been developed over the past few decades.
Patterson and Hamblin (1988) developed a snow and one-component ice version (DYRESMI) of the 1D Dynamic Reservoir Simulation Model (DYRESM; Imberger and Patterson, 1981) that couples to a thermodynamic lake-mixing model of the water column beneath the ice. Gosink (1987) independently developed an ice-cover routine for DYRESM (DYRESMICE). These formulations of ice cover were suitable for high-latitude lakes, but lacked important processes necessary for mid-latitudes. Gu and Stefan (1993) proposed a two-component ice and snow model that included sediment heat release and the melting of ice and snow by rain. However, the snow and ice parameters were considered constant. This model was later modified by Fang and Stefan (1996) for sediment heat release and vertical thermal diffusivity to better suit deeper lake basins.
The Mixed Lake with Ice (MLI) model of Rogers et al. (1995) extended DYRESMI to include snowmelt by rain, the formation of white ice, and the variability of snow density and albedo in order to address particular concerns at mid-latitudes. However, the underlying water-column dynamics were not coupled to the ice model in their study, because it was developed for shallow lakes with negligible vertical thermal structure. The bulk water temperature was simply predicted as a function of the solar radiation reaching the top of the water column. Similarly, McCord and Schladow (2000) developed the Dynamic Lake Model (DLM) for the study of artificial aeration kinetics in ice-covered lakes. The hydrodynamic component of DLM is similar to DYRESM and the ice cover sub-model is based on MLI.
More recently, a pseudo two-dimensional version of DYRESM was developed by Yeates and Imberger (2003) with explicit benthic boundary layer and internal mixing algorithms. The objective of the present study is to couple the 2003 version of DYRESM with a new three-component ice algorithm which is similar to MLI but employs two-way coupling between the ice/snow layers and the underlying water column. This allows year-round simulation of lake and ice dynamics in deep lakes at temperate latitudes. The results are validated against published observations from Harmon Lake (British Columbia) and new observations from Eagle Lake (Ontario). These comparisons show the model to successfully simulate snow and ice thickness, winter stratification, and spring turnover. This model has the advantage of being coupled with the advanced biogeochemical Computational Aquatic Ecosystem Dynamics Model (CAEDYM; e.g., Antenucci et al., 2003; Trolle et al., 2008) and will be used in future to determine the potential impacts of climate change on lake biogeochemistry.
METHODS
Model description
DYRESM is a 1D vertical Lagrangian layered model. The Lagrangian layers dynamically adjust their thicknesses according to the development of stratification, turbulent mixing and/or changes in lake volume. This eliminates numerical diffusion, which acts to artificially thicken the thermocline in Eulerian models (Laval et al., 2003). DYRESM is based on the assumption that density stratification slows vertical mixing, causing vertical scalar gradients to persist, while horizontal gradients are quickly mixed by diffusion, advection, and convection. Technical details on DYRESM can be found in Imerito (2007) and Yeates and Imberger (2003). In the present study, an ice formation algorithm was implemented in the DYRESM surface layer thermodynamics routine. This algorithm solves the quasi-steady heat transfer equation in three layers (blue ice, white ice, and snow) located between the water and the air (Fig. 1).
The processes included in this algorithm include the penetration of shortwave radiation through the snow and ice layers to the water column below. Non-penetrative longwave radiation contributes only to the melting of snow or ice and does not reach the water surface. Rain-induced melting, the formation of blue and white ice, and the variability of snow density and snow and ice albedo are also included in the ice formation algorithm, as these parameters vary with mid-latitude weather patterns. Blue ice initially forms when temperatures are below freezing at the water surface. Snow then accumulates on top of the blue ice. The accumulated snow floods when the blue ice cannot support the snow weight, and subsequently white ice forms.
These processes were parameterized into the steady state equation for heat transfer, which was solved for the contiguous layers of blue ice, white ice and snow, between the atmosphere and water column (Fig. 1). We neglect advection in the heat conduction equation as the Stefan number S=ci T/Lw <0.1, where Lw is the latent heat of fusion of water, ΔT is the temperature difference across the ice, and ci is the specific heat capacity of ice. In other words, the temperature distribution reaches steady state within one time step (Gilpin et al., 1980; Hill and Kucera, 1983). Therefore, the steady state heat conduction equation can be written as (a list of notations is provided in Tab. 1):
(eq. 1)
where T is the temperature, Q the heat flux, KH the diffusivity of heat, ρn the snow or ice density, cn the specific heat of snow or ice, and z the vertical coordinate. Combining this with the penetrating radiation I=Isw exp(-λ(h-z)), thermal conductivity K=Khpncn over the three distinct layers (i = blue ice, e=white ice and s=snow) in Fig. 1 gives
(eq. 2)
With the following boundary conditions
(eq. 3)
Here, Isw is incoming solar radiation equal to 1-αQsw, where α is the albedo of the ice and snow cover. Isw has two components (Boer, 1980; Kirk, 1983): visible (subscript 1, A1=0.7) and near infra-red (subscript 2, A2=0.3), with a distinct attenuation coefficient λ for each component and each ice and snow layer. Qsi is the heat flux per unit volume released during the formation of white ice that results from the flooding of the snow layer. Tf and Ta denote the freezing and air temperatures, respectively. Eq. (2) with the eq. (3) boundary condition can be written as:
(eq. 4)
With the following value of the heat flux at the interface between air and ice or air and snow cover, eq. (4) can be solved
(eq. 5)
In eq. (4), Q0, Qlw, Qs, Qe, and Qr are the incident shortwave solar radiation, the net longwave radiation, the sensible heat flux, the latent heat flux, and the heat flux due to rainfall, respectively. If Qlw is not provided as model input, it can be calculated within DYRESM from cloud cover observations (Swinbank, 1963). Qs and Qe are calculated as in the absence of ice (Fischer et al., 1979) with the corresponding snow/ice cover temperature instead of the water surface temperature and the latent heat of vaporization for ice.
Lw, Qp, Qf, and Qw are the latent heat of fusion of water, the penetrative shortwave solar radiation, the heat flux from the ice to the water and the heat flux from the water to the ice, respectively. Ice production and melting is the result of an imbalance between Qf and Qw (dhi/dt=(Qf-Qw/piLw).
The ice formation algorithm was implemented in the DYRESM surface layer thermodynamics routine. To simulate the temperature and ice cover during winter in midlatitude lakes, with relatively rapid changes in weather conditions, other important processes such as variable density and conductivity of snow, variable snow and ice albedo, formation of white ice and melting of ice and snow by rain have been included in the DYRESM ice formation model. In the following, we describe these important processes.
Snow density
The density of newly fallen snow depends on the air temperature and on changes resulting from compaction. Weather can change significantly in mid-latitude regions over the time-scales associated with synoptic weather patterns. Therefore, variablity in process-based parameters, such as snow density, are important. Conversely, the density of newly fallen snow in high-latitude regions, where the winter temperature is consistently cold, is nearly constant. Snow compaction can also cause significant variations in snow density with time. We consider the change in snow density based on Koren et al. (1999):
(eq. 6)
where B=tC1exp(C3Ts-C2ρs), ρsc is the compacted snow density, ws=hsρs/ρw is the water equivalent of the snow, C1=2.77x10−4m1s1, C=0.021 m−3 kg and Q=0.08°C−1.
The density of freshly fallen snow is a function of air temperature (Gottlieb, 1980), whereas in MLI freshly fallen snow has a constant density of 80 kg m−3. The density of fresh snow is given by ρs,new = max[50,1.7(Ta+15)1.5], implying that the density approaches a minimum of 50 kg m−3 as Ta approaches -5°C. The average snow density is calculated from
(eq. 7)
where ws,new and ws,c refer to the relative water equivalent of the newly fallen snow and of the compacted existing snow, respectively.
Snow conductivity
The thermal conductivity of snow increases with increasing snow density (Ashton, 1986):
(eq. 8)
Snow and ice albedo
Recent numerical modeling of lake ice (Yang et al., 2012) has confirmed the importance of snow and ice albedo. In mid-latitude lakes, where the weather often changes rapidly, it is important to consider the variation of snow and ice albedo, as this regulates the penetration of shortwave solar radiation, which is an important source of energy driving the thermal structure in insulated ice-covered lakes (Kirillin et al., 2012). The albedo of snow and ice (αs and αi, respectively), were calculated following Vavrus et al. (1996). The albedo decreases with the thickness of blue ice and snow (hi and hs, respectively), as well as with increases in ice or snow temperature up to the melting point. For example, the snow albedo decreases with the thickness of the snow to account for non-even distributions of ice and snow, and consequently the areas without snow and ice cover. The albedo equations are
(eq. 9)
(eq. 10)
The model first calculates the albedo of ice, regardless of the existence of a snow layer, because αi is required in the calculation of snow albedo when the snow thickness is less than 0.1 m. The snow albedo is then calculated if required. Here, we have assumed that both blue ice and white ice have the same albedo due to the difficulty in determining the albedo of white ice, which is typically higher than that of blue ice.
White ice
In the model, white ice forms when the blue ice is unable to support the weight of accumulated snow, which becomes submerged and subsequently freezes. This occurs if hs>hs,max, where hs,max is the maximum snow thickness that can be supported by the existing blue ice
(eq. 11)
The increase in the white ice thickness (Δhe) and in the heat released (Qsi) are (Rogers et al. , 1995):
(eq. 12)
and
(eq. 13)
Then hs=hs,max and he=he+he. This parameterization suggests immediate flooding and instant formation of white ice when eq. 11 is satisfied. However, the blue ice may hold the weight of snow at hs>hs, max for days or weeks before cracks form and water flows on to the ice surface. The water can also remain unfrozen under the snow, without white ice formation, for up to several weeks if the snow cover is thick and/or air temperatures are high. These limitations can only be minimized by applying lake-ice models and adjusting the governing equations given the limited availability of observed data.
Rain
(eq. 14)
where R is the amount and ca the heat capacity of rain.
The other fixed model parameters used in this study are presented in Tab. 2.
Validation of the ice algorithm
The modified version of DYRESM was applied during the winter season to simulate the water column temperature profile and ice and snow accumulation in Harmon Lake, British Columbia, and Eagle Lake, Ontario (Fig. 2; Tab. 3). These lakes have simiar depths, but are in two distinct mid-latitude regions of Canada and have significantly different bathymetric shapes. Eagle Lake has a surface area 26 times greater than Harmon Lake and has a complex topography with many islands. Moreover, the lakes are in different climactic regions of Canada. Eagle Lake is in central Canada on the Canadian Shield, whereas Harmon Lake is an alpine lake in western Canada. Given these differences, we feel the model has been suitably validated. The Eagle Lake simulation was designed to test the ability of the model to capture a whole year of simulation, including ice-on and ice-off periods, while the Harmon Lake simulation allows the modeling of the accretion and ablation of ice during a period of ice cover to be tested.
Application to Harmon Lake
A comprehensive field study was conducted by Rogers et al. (1995) from 13 December 1991 to 24 March 1992. The temperature was measured at intervals of 1 m by profiling down to the maximum depth of the lake. Snow and ice thickness were measured biweekly, with no distinction made between blue and white ice. Meteorological data (Fig. 3), including shortwave solar radiation, wind speed, relative humidity, and air temperature, were measured at Menzies Lake, about 5 km northeast of Harmon Lake and 100 m higher in elevation. Longwave radiation was estimated from cloud cover (Swinbank, 1963) within DYRESM. Cloud cover observations (as reported in Rogers, 1992) were determined from the difference between the observed and theoretical clear-sky shortwave radiation (Tennessee Valley Authority 1972). Rainfall was measured at Merritt, 15 km northwest of Harmon Lake, and snowfall was measured at Menzies Lake using a snowboard with an accuracy of ±1 cm after each snow event; however, reported errors are much greater (Rogers et al., 1995) due to significant snow compaction. Because of potential inaccuracies in precipitation measurements related to local conditions (e.g., snow drift), and the strong influence of precipitation on ice formation, observed snow and rainfall values at Menzies Lake were compared with the observed snow thickness at Harmon Lake, and the adjusted values of Oveisy et al. (2012) were used to compensate for snow drift and other inconsistencies. In small lakes, sediment heat flux is typically 1-5 Wm−2; here we follow the suggestion of Rogers et al. (1995) and use a constant 2.8 Wm−2 during the simulation. The model was set up and run from 13 December 1991 to 24 March 1992. The initial ice and snow thicknesses and temperature profile were taken from the observations on 13 December.
Harmon Lake model results
Fig. 4 shows the temporal evolution of the observed and modeled snow and ice thicknesses. The initial modeled and observed snow depths decrease over the first 30 days due to compaction, then increase rapidly due to heavy snowfall (days 25-38, Fig. 3b) and then melt during heavy rain (days 46-50, Fig. 3b). The increase in ice thickness during winter and subsequent melting during spring are simulated well. The ice thickness increases from 0.15 m to 0.35 m during the first 80 days, with the rate of ice growth being sensitive to the depth of overlying snow. For example, the snow on days 40-46 insulates the ice surface, causing a decrease in the rate of ice growth to near zero, and a subsequent increase in ice growth (days 51-75) after the rain-induced snowmelt near day 51. The ablation of ice from days 75-85 is slow, and with an increase in solar radiation and air temperature the ice starts to melt rapidly after day 85. The model correctly responds to this rapidly changing meteorological forcing, which is characteristic to mid-latitude regions.
The water column beneath the ice responds to the ice-induced changes in surface thermodynamic fluxes (Fig. 5). The modeled volume-averaged lake temperature shows a slight decrease in the first 15 days (Fig. 6) as the total heat flux from water to ice is greater than the sum of the solar radiation reaching to water and the sediment heat flux, 2.8 W m−2 (Fig. 5). Thereafter, the lake temperature increases when the sediment heat flux is trapped due to the insulating effects ofthe freshly fallen snow (day 40-70, Fig. 6) and the decrease in the water-to-ice heat flux (Fig. 5). Snowmelt and increases in air temperature and solar radiation cause the average lake temperature to increase rapidly during the last 40 days of the simulation. The disappearance of the snow allows more solar radiation to penetrate the ice (ice albedo < snow albedo), increasing the temperature below. To a lesser degree, the increased air temperature increases the temperature gradient and the resulting temperature flux through the ice layer. Profiles of modeled and observed temperatures (Fig. 7) show the characteristic inverse winter stratification to persist until day 90, when the spring turnover mixes the water column in response to the increased air temperature and solar radiation. The sharp temperature gradient occurs under the ice cover from 3°C to 0°C in both the model and the observations. The comparison of modeled with observed temperature profiles is excellent, with the Mean Bias Deviation (MBD)=-2.8% and the Root Mean Square Deviation (RMSD)=0.38°C, although the surface water in the model warms later than is observed. This likely results from the laterally averaged model not being able to reproduce differential heating in shallow, nearshore regions relative to the deeper, mid-basin site.
The main heat loss from the lake is through the longwave heat flux (Fig. 5), while the sensible heat flux is almost always into the lake, as the ice or snow temperature calculated by solving the heat conduction equation is always lower than the air temperature (Rogers et al., 1995). The latent heat flux is negligible during the entire simulation (>0.7 W m2). During the first 60 days, the heat flux from the ice to the water is <10 W m2, particularly during days 12-50. In this period the heat flux from sediment is important, as its magnitude (2.8 W m2, not shown in Fig. 5) is similar to the heat flux from the water to the ice, while the penetrative shortwave solar radiation reaching the water is also <5 W m2. Comparing the penetrative solar radiation reaching the water with the incoming solar radiation (Fig. 2) reveals that during thick ice and snow cover only small portion (~20%) of the solar radiation reaches the water, and this increases thereafter to ~60% as the snow thickness decreases toward the end of the simulation. The heat flux from the water to the ice increases near the end of the simulation because of the increased water temperature beneath the ice that results from the increased penetrative shortwave radiation. The rain-induced heat flux is significant during heavy rainfall (days 46-50, Fig. 3b) and the sensible and latent heat released from this rainfall decrease the snow depth significantly.
Comparisons to MLI and DYRESMI (from Rogers et al., 1995) are shown in Figs. 4 and 6. In DYRESMI the initial predictions of ice and snow thicknesses are not the same as observed (Fig. 4), because the initial ice thickness could not support the weight of the observed snow due to the high snow density (constant at 330 kg m−3) assumed in DYRESMI (Rogers et al., 1995). DYRESMI also overestimates the thickness of the ice by ~50% relative to the observations during most of the simulation period (Fig. 4). This is due to the high snow conductivity in DYRESMI (constant at 0.31 Wm *°C−1). The snow depth is also predicted to be higher than observed, likely because of the assumption of no settling and no snow drift, and because snow melt due to rain is neglected. MLI underestimates ice thickness by ~11% and snow thickness by ~10% compared to the observations, because the water temperature is approximated from the solar radiation reaching the water, and the water column dynamics beneath the ice are neglected. The prediction of lake-average temperature using DYRESMI is significantly lower than observed (Fig. 6). This is due to an overestimation of the snow and ice thickness, which results in excessive insulation of the lake and neglects the heat flux released from the sediment. MLI predicts the average lake temperature accurately for the first 50 days, but overestimates the temperature at which the snow melts on the ice, because ice has a much smaller albedo (0.02 in Rogers et al., 1995) than snow (~0.80 in Rogers et al., 1995), causing solar radiation to warm the lake quickly. The albedo of ice is constant in MLI, while it is variable in our formulation.
To further investigate the effects of using constant snow and ice albedos for the simulation of temperate lake winter dynamics, we ran the model with constant albedo values (Fig. 4). For thick ice and snow with cold air temperatures (<-5°C) and high albedo (ice and snow albedos of 0.6 and 0.7, respectively). Fig. 4 shows that these values result in an overestimate of ice thickness resulting from the excessive insulation of the lake from solar radiation. The overestimate increases toward the end of the simulation, when the model-variable albedo decreases significantly with the decrease in simulated ice and snow thickness and the increase in ice and snow temperature. For mid-rage constant albedos (an ice albedo of 0.44 and a snow albedo of 0.50), which are characteristic of thick ice and snow with the temperature close to melting point, the model predicts the ice thickness well during relatively thick ice and snow conditions (before day 53), while overestimating it thereafter, as the model does not capture the decrease in the ice/snow albedo with decreasing snow/ice thicknesses. For runs with low ice and snow albedos, characteristic of the thin ice and snow layers (0.1 m and 0.05 m ice and snow thickness with an ice albedo of 0.31 and a snow albedo of 0.40, respectively), the ice thickness is underestimated for the entire simulation, but the rate of ice melting is similar to that with variable ice/snow albedo. This implies that the values are valid near the end of the simulation, when the ice and snow thicknesses are small.
Application to Eagle Lake
Eagle Lake (Fig. 2b, Tab. 3) has a complex bathymetry and is located 65 km north of Kingston, Ontario. Although oligotrophic, the lake has a thin hypolimnion, and so oxygen is rapidly consumed by the natural sediment oxygen demand (Bouffard et al., 2013). Consequently, lake trout are threatened and there is a moratorium on shoreline development.
The water-column temperature was measured at St. 1 and 2 in Eagle Lake (Fig. 2b) using PME T-chains (www.pme.com) with sensors (accuracy of 0.010°C) sampling at 0.1 Hz every 0.5 m from 2 m below the surface to the lake bed. The logger at St. 1 failed soon after deployment. Ice and snow thickness were measured manually during two winter field trips, with significant local variation in snow thickness observed to result from wind drift. To force the model, wind speed and direction, incident shortwave solar radiation, relative humidity, and air temperature were measured on an exposed island north of St. 2 using a HOBO U30 weather station mounted on a 3m tripod. The three-day average solar radiation, air temperature and precipitation are presented in Fig. 8. As in Harmon Lake, cloud cover was computed from the difference between measured and clear-sky shortwave radiation (Reed, 1977) and longwave radiation was calculated from cloud cover within DYRESM (Swinbank, 1963). Precipitation was obtained from the Environment Canada Hartington Station, which is ~25 km south of Eagle Lake. The model was started on 16 May 2011 (day of the year: 137) with initial conditions from the in-lake measurements, and run for 374 days until 25 May 2012 (day of the year: 146). As the precipitation observations were not on site, and to account for snow drift, we used the same approach as for Harmon Lake and adjusted the snowfall to match the on-site measurements.
Eagle Lake model results
The model shows good skill in predicting ice and snow in comparison to the limited manual measurements during the two field trips (Fig. 9). Modeled ice forms on day 211 (14 December 2011) with a rapid increase in the ice thickness thereafter until day 224 (28 December, 2011). As a result of several snow events between days 224 and 265 (Fig. 8b), snow cover forms thereafter (Fig. 9) and the rate of the ice growth decreases as the ice becomes insulated. When the snow cover melts (days 253 and 265), a sharp increase in ice thickness occurs. The ice thickness reaches a maximum near 45 cm on day 277 (19 February 2012). The ice then slowly melts as air temperatures increase into spring, with the melt-rate increasing after the snow cover melts. After snowmelt, the ice thickness sharply decreases and disappears by day 325 (7 April 2012). Observations of ice-on and ice-off from nearby Little Silver Lake (15 km NW of Eagle Lake) indicate that ice typically forms in late November to late December and disappears by early April (Boegman et al., 2012), which is in agreement with the model results that predict ice-on to have occurred on 14 December 2011, and ice-off on 7 April 2012.
The modeled temperature profile follows the observations well throughout the year-long simulation (Fig. 10). Both modeled and observed temperatures increase to ~25°C by mid-summer on day 69 (26 July 2011) and then begin to decrease on day 85 (11 August 2011). Fall and spring turnover are recognizable around day 205 (9 December 2011) and day 333 (14 April 2012), respectively. Due to the absence of wind stress during the period of ice cover, the characteristic inverse winter stratification and sharp temperature gradient immediately beneath the ice are both modeled. To prevent the mooring freezing into the ice, no observations were made in this region of the lake, which is within 2 m of the surface (Fig. 10b). Due to the insulation of the lake by the ice and snow layers, the temperature profile does not significantly vary during the period of ice cover. The modeled and observed water column warms quickly after day 335 (17 April, 2012) due to melting of snow and the lake surface temperature reaches 15°C by the end of the simulation on day 374 (26 May 2012).
To investigate the importance of variable snow and ice albedo for Eagle Lake, we also ran this model application with constant values (Fig. 9). The results are similar to those from Harmon Lake. With high albedos (0.7 and 0.6 for snow and ice, respectively), the ice thickness was overestimated and the lake was ice-covered until day 354 (6 May 2012). Using lower constant values (0.50 for snow and 0.44 for ice), the prediction of ice thickness is similar to that with the original variable ice and snow albedo routines, until day 235 (8 January, 2012), when the thickness of ice and snow spuriously increases when air temperatures are below -5°C (Fig. 8a); subsequently, this run underestimates the ice thickness. As the air warms toward spring and snow disappears, the constant albedo causes a higher rate of increase of ice thickness, relative to the variable routine, and after day 305 the model overestimates ice thickness. For low albedos of 0.31 for ice and 0.405 for snow, respectively, the model underestimates the ice thickness from day 220.
DISCUSSION
The process interactions of the water column and ice algorithm are captured in the current model with reasonable accuracy. In small lakes, such as Harmon Lake, this interaction is more recognizable as these lakes quickly react to external heat fluxes. We have shown that for temperate mid-latitude lakes, capturing these interactions requires modeling the effects of transient mid-latitude meterological processes, such as snowmelt due to rain and variation of snow and ice albedo, because these lead to relatively rapid changes in the characteristics of the snow and ice layers. Changes in these layers ultimately influence the heat fluxes to the underlying water column and the resultant winter circulation dynamics (Wang et al., 2010; Oveisy et al., 2012; Yang et al., 2012; Kirillin et al., 2012); therefore, accurate estimation of these parameters is very important. The model accuracy in Harmon Lake (ice thickness, MBD=-1.1% and RMSD=0.01 m; snow thickness, MBD=-3.1 % and RMSD=0.02 m) is higher than in previous 1D simulations of lake ice (e.g., Rogers et al., 1995: ice thickness, MBD=-10.6% and RMSD=0.06 m; snow thickness, MBD=-9.6% and RMSD=0.03 m; Patterson and Hamblin, 1988: ice thickness, MBD=52.9% and RMSD=0.14 m; snow thickness, MBD=22.3% and RMSD=0.07 m). This is because of the more accurate process-based ice formulation algorithm in the present model (i.e., variable ice and snow albedo, variable snow density and conductivity), the inclusion of three snow and ice layers, and two-way coupling between the ice and the water column. The heat exchange between Harmon Lake and the air is also predicted more accurately in this model, relative to previous models, as evidenced by the better agreement between predicted and observed mean lake temperature (present model: MBD=-1.1% and RMSD=0.1°C; Rogers et al., 1995: MBD=5.6% and RMSD=0.3°C; Patterson and Hamblin, 1988: MBD= -13.4% and RMSD=0.57°C).
The current model was also successful in simulating year-round temperatures in Eagle Lake (RMSD=1.8°C and MBD=8.80%). The lower accuracy compared to Harmon Lake was expected due to the difficulty in representing the complex bathymetry of Eagle Lake with a 1D model (Fig. 2) and the Harmon Lake model being run for a significantly shorter time period, during which there was complete ice cover. A similar reduction in accuracy was also observed by Fang and Stefan (1996) in year-round modeling, compared to that during the period of ice-cover.
The present model shows a similar skill to other models (Fang and Stefan, 1996) in reproducing observations of temperature profiles (RMSD=1.1-1.3°C) and ice thicknesses (RMSD=0.07-1.2 m) in Thrush Lake, Minnesota, and Little Rock, Wisconsin, emphasizing that model errors may not be dependent so much on variations in general model structure (e.g., mixing routines) but rather on inaccuracies in estimating meteorological forcing and/or the mathematical definitions of particular snow and ice parameters (e.g., albedo). Three-dimensional modeling of Harmon Lake using ELCOM (Oveisy et al., 2012) showed similar accuracy in predicting ice thickness (RMSD=0.01 m), snow thickness (RMSD=0.02 m), and temperature profiles (MBD=1.23% and RMSD=0.40°C) to the current 1D model using similar ice formation algorithms, confirming that in small lakes with simple geometries, a 1D vertical model can estimate ice and snow thickness and, therefore, inter-annual temperature dynamics, as accurately as a 3D model. This has implications for the inclusion of the many small mid-latitude lakes in climate models (Mackay, 2012).
CONCLUSIONS
An ice formation algorithm, consisting of three layers of blue ice, white ice and snow, was added to DYRESM. The model was designed to capture the variation in snow and ice formation processes occurring at mid-latitudes, where rapid changes in weather conditions occur. These processes included snow and ice melting due to rain, variable ice and snow albedos, and variable snow conductivity. With the implementation of this algorithm, DYRESM can now be run for the multi-year and long-term prediction of lake dynamics (e.g., to investigate lake response to climate change). The model was validated against observations from two lakes with an accuracy of ~1 cm in the prediction of ice thickness and <2°C in the prediction of water temperature, and showed superior results compared to previous models due to the inclusion of the effects of meterological variability and two-way coupling between the overlying ice and the water column. The results show the ice formation algorithm to be highly sensitive to the snow and ice parameters. Further winter observations of snow and ice properties are needed to improve model accuracy and applicability. These simulations will eventually be extended to include impacts of climate on lake biogeochemistry following the adaption of CAEDYM to winter conditions.
Abstract
Main Text
INTRODUCTION
METHODS
Model description
Snow density
Snow conductivity
Snow and ice albedo
White ice
Rain
Validation of the ice algorithm
Application to Harmon Lake
Harmon Lake model results
Application to Eagle Lake
Eagle Lake model results
DISCUSSION
CONCLUSIONS