2. Basic theory and development of a simple model

The growth of an ice cover on a body of water is essentially described by the Stefan Problem. Based on the original description of the physical process in mathematical terms (Stefan, Reference Stefan1891), many studies have been published on analytical and numerical solutions of the Stefan Problem (e.g. Foss and Fan, Reference Foss and Fan1972; Hutter and others, Reference Hutter, Zryd and Röthlisberger1990). To provide a background, we briefly review the simplest possible model for ice growth on a lake (Ashton, Reference Ashton1989). It is assumed that the ice surface temperature T _s equals the air temperature, and that there is no snow. The assumption of no snow cover, in particular, is unrealistic for alpine lakes, but the case may still serve as a reference case which provides an upper bound of ice growth.

For a top temperature T _s and bottom temperature T _m (melting point, set to 273.15 K), the diffusive heat flux Q through the ice layer in steady-state conditions is given by

(1)$$Q = k_{\rm i}( {T_{\rm s}-T_{\rm m}} ) /h_{\rm i}.$$

here h _i is the ice thickness and k _i is the thermal conductivity of the ice. At the bottom of the ice layer the heat flux is balanced by the heat of fusion released during the freezing process. If the growth rate of the ice layer is denoted by dh _i/dt, we have

(2)$${\rm d}h_{\rm i}/{\rm d}t = \displaystyle{{k_{\rm i}} \over {\rho _{{\rm i\;}}L}}( {T_{\rm m}-T_{\rm s}} ) /h_{\rm i}, \;$$

where ρ _i is ice density (917 kg m³) and L is the latent heat of fusion ($334\;{\rm kJ}\;{\rm k}{\rm g}^{ \hbox{-} 1}$).

For a constant negative temperature difference (T _s − T _m), this can be integrated to give

(3)$$h_{\rm i}( t ) = \sqrt {h_{{\rm i}, 0}^2 + \displaystyle{{2k_{\rm i}} \over {\rho _{{\rm i}\;}L}}( {T_{\rm m}-T_{\rm s}} ) \;( {t-t_0} ) {\rm \;\;}} , \;$$

where h _i,0 is the ice thickness at initial time t ₀. This solution has been known for a long time, and it is quite fundamental because it shows that ice growth is basically parabolic in time. This is due to the decreasing temperature gradient and smaller associated heat flux in the ice when the ice layer thickens.

The analysis above shows that the surface temperature of the ice layer is a key variable that determines the rate of thickening. T _s can be calculated either as a skin temperature or as a temperature of a finite upper layer with constant or varying heat capacity. In advanced lake ice models all energy fluxes are normally considered on a vertical grid of variable resolution, in which new ice or snow layers may be generated or removed. In our approach we do not deal with vertical resolution. Instead we include a simple scheme based on the assumption of a quasi steady-state heat flow of equal magnitude in the ice and snow layer, yielding an expression for T _s in terms of air temperature, snow depth and the thicknesses of the ice and snow layer. Our scheme implies a huge simplification, but the resulting model is capable of predicting ice growth rates with good accuracy. Before formulating the model, we list a few considerations that are at the basis of our approach:

• • During the early/midwinter growth of ice on Lake St. Moritz, air temperature is the most important factor. Solar elevation is low, implying a small energy input, especially once the ice is covered by snow.

• • It is well known that clear nights facilitate ice growth because the longwave radiation balance is strongly negative. However, clear nights over a wide horizontal surface have very low associated air temperatures, implying that the effect of longwave cooling can be dealt with to a large extent through the air temperature.

• • Daily and short-term temperature variations damp out rapidly with depth in the ice layer (as demonstrated by for instance Gould and Jeffries (Reference Gould and Jeffries2005), who performed detailed temperature measurements during lake ice formation). We therefore construct a model version with daily input.

• • As has been demonstrated in several studies, snow is the most obvious factor affecting ice growth (e.g. Fichefet and Maqueda, Reference Fichefet and Maqueda1999; Yang and others, Reference Yang, Leppäranta, Cheng and Li2012). The inclusion of the insulating effect of snow is essential even in the simplest model.

• • Mechanical grooming of snow has been shown to be effective in cooling the underlying ice or soil (e.g. Keller and others, Reference Keller2004). It is also applied over a large area (about 1/3) of Lake St. Moritz, which facilitates growth of the ice layer.

We now turn to the description of Simple Lake Ice Model (SLIM). We calculate T _s according to

(4)$$\displaystyle{{{\rm d}T_{\rm s}} \over {{\rm d}t}} = ( {T_{\rm s}^\ast{-}T_{\rm s}} ) /\tau , \;$$

where $T_{\rm s}^\ast$ is an equilibrium ice surface temperature determined by the air temperature and the snowpack. Eqn (4) describes how the ice surface temperature evolves to $T_{\rm s}^\ast$, with an exponential response timescale τ (typically a few days). To find an expression for $T_{\rm s}^\ast$ we assume that the temperature at the ice or snow surface is always equal to the air temperature T _a, and that the diffusive heat flux through the ice layer and through the snow pack are equal. This implies

(5)$$\displaystyle{{k_{\rm i}( {T_{\rm m}-T_{\rm s}^\ast } ) } \over {h_{\rm i}}} = \displaystyle{{k_{\rm s}( {T_{\rm s}^\ast{-}T_{\rm a}} ) } \over {h_{\rm s}}}, \;$$

where h _s is the snow depth and k _s is the characteristic thermal conductivity of the snowpack (which may depend on the mean density). Solving for the ice surface temperature yields

(6)$$T_{\rm s}^\ast{ = } \displaystyle{{h_{\rm s}T_{\rm m} + rh_{\rm i}T_{\rm a}} \over {h_{\rm s} + rh_{\rm i}}}, \;$$

where

$$r = k_{\rm s}/k_{\rm i}.$$

When there is no snow we thus have $T_{\rm s}^\ast{ = } T_{\rm a}$. When the snow depth becomes very large, we have $T_{\rm s}^\ast \to T_{\rm m}$, implying that the ice layer becomes isothermal and the heat flux approaches zero. The ratio of the thermal diffusivity of the snowpack to that of ice will always be larger than 1, and can easily have a value of 10 or even 20 for fresh dry snow (e.g. Sturm and others, Reference Sturm, Holmgren, König and Morris1997).

The prognostic equation for the ice thickness is similar to Eqn (2), but slightly modified by introducing a parameter δ :

(7)$${\rm d}h_{\rm i}/{\rm d}t = \displaystyle{{k_{\rm i}} \over {\rho _{{\rm i\;}}L}}( {T_{\rm m}-T_{\rm s}} ) /( {h_{\rm i} + \delta } ) .$$

The introduction of δ prevents the occurrence of a very high heat flux when the ice thickness is very small (which could cause instability of the numerical scheme). Including δ also provides some control on the heat flux through the ice, which may deviate from the value of pure ice due to the presence of impurities and air bubbles. In addition to this, the freezing of the lake is not a uniform process. Some parts may freeze one or two days earlier than other parts, related to differences in exposure and wind conditions. We are not able to model this initial stage of ice formation, but will use the parameter δ to get a good match between simulated and observed ice thickness as soon as the first measurement can be done.

The model formulated above has three parameters (r, τ and δ) that can be used as tuning parameters. The numerical implementation of the model is very simple. Eqn (7) is easily integrated with a forward Euler scheme with a time step of 1 h. The meteorological forcing is through the prescription of the daily mean air temperature, which is entered into Eqn (5). As initial condition a certain ice thickness different from zero should be given (e.g. h _i = 0.02 m) at a certain day.

3. Application to the winter of 2021/22

It is common practice to measure ice thickness at many sites on the lakes of the Upper Engadin. However, as soon as the thickness reaches 0.3 to 0.4 m, which normally occurs in the course of January, the measurements are stopped because the ice is considered safe enough. Fortunately, in the winter of 2021/22 the observations on a central part of Lake St. Moritz were continued until mid-February. The lake froze over entirely on 16 December, and ice thickness and snow depth measurements were carried out at fixed locations on 19 December, 27 December, 2 January, 15 January, 28 January and 13 February. At each location three to four measurements were done within a distance of about 10 m to obtain a representative value. The variability among these measurements was quite small (typically 5%).

Daily mean air temperature (T _a in Eqn (6)) has been taken from the SwissMeteo weather station at Samedan airport, about 4 km away from the lake. Although the altitude of the weather station is about 60 m lower, we have not made any corrections because there are no indications that wintertime conditions in St. Moritz are colder than in Samedan. Although it may be windy occasionally, rather calm conditions are the rule and the snow cover over the lake is normally quite uniform. Rain events and blizzards are very rare. At the sites of the ice thickness measurements the snow was always groomed by machines to facilitate transportation of material and to increase ice growth. Unfortunately, snow density measurements have not been carried out, but it is known that density of groomed snow has a typical value of $500\;{\rm kg}\;{\rm m}^{ \hbox{-} 3}$ (Keller and others, Reference Keller2004).

For the parameter τ we have chosen a value of 2.5 days. After some numerical experimentation this value appeared to broadly reproduce the ice temperature measurements as reported in the study by Gould and Jeffries (Reference Gould and Jeffries2005). The parameter δ was determined by matching the first ice thickness measurement with the calculated ice thickness, yielding δ = 0.09 m. During the initial freezing period there was no snow on the ice.

The optimal value of the parameter r, which determines the insulating effect of the snow cover, was found by minimizing the root-mean-square difference σ between simulated and observed ice thickness, yielding r = 4.9 and σ = 2.1 cm.

The result of the simulation is shown in Figure 2. The match between observed and simulated ice thickness is quite good, although the model slightly overestimates the slow-down of the growth rate after the first snowfall. The delayed response of the ice surface temperature to air temperature variations and snow cover is clearly seen. From this simulation we conclude that the essential processes are captured by SLIM. However, we realize that the conditions were relatively favourable: not too much snow and no really warm spells with air temperatures well above the freezing point for several days.

Fig. 2. Simulation of ice thickness on Lake St. Moritz, for the period 1 December 2021–28 February 2022. Ice surface temperature and air temperature are shown in red (scale at left). Daily values of calculated ice thickness are shown in blue (scale at right). Observed ice thickness is shown by black diamonds (dates given by the labels). The green curve shows the observed snow depth on the ice, which serves as input for the model calculation. On 21.2 an additional snow depth measurement was made.

We have also applied the model to Lake Silvaplana, which is about 5 km southwest of Lake St. Moritz. Unfortunately, ice thickness and snow depth measurements have been made only on 27 December, and 3, 13 and 17 January (ice thickness 34 cm). The lake is twice as deep as Lake St. Moritz, which probably is the reason that it froze over on 24 December (as compared to 16 December for Lake St. Moritz). Without changing the values of the model parameters, the model simulation was reasonably good (σ = 3.8 cm), although the simulated growth after day 40 is slightly larger than observed. The simulation can be improved by decreasing the thermal conductivity of the snow. For r = 15 the simulated curve is in between the observations from 13 and 17 January (also shown in Fig. 3), and we have σ = 2.3 cm.

Fig. 3. Simulation of ice thickness on Lake Silvaplana, for the period 1 December 2021–28 February 2022. Ice surface temperature and air temperature are shown in red (scale at left). Daily values of calculated ice thickness are shown in blue (scale at right). Observed ice thickness is shown by black diamonds (dates given by the labels). The green curve shows the observed snow depth on the ice, which serves as input for the model calculation. The curve labelled ‘large r’ shows ice thickness for the simulation with the strongly reduced thermal diffusivity of the snow.

We have done some sensitivity tests to investigate more explicitly the role of snow in the model. Two extreme scenarios were simulated: (i) no snow cover at all, and (ii) snow that is not groomed and keeps a density of 200 kg m⁻³ throughout the period of simulation. Measured thermal conductivity of snow varies enormously, depending mainly on the structural properties of the snow. Following Sturm and others (Reference Sturm, Holmgren, König and Morris1997), for a density of 200 kg m⁻³ we have used as a characteristic value of $k_{\rm s} = 0.1\;{\rm W}\;{\rm m}^{ \hbox{-} 1}\;{\rm K}^{ \hbox{-} 1}$, corresponding to r = 22.

The result is shown in Figure 4. At day 90, the ice thickness for the case without snow would be 0.73 m, with low density snow it would be only 0.45 m. Grooming of the snow clearly helps to makes the ice thicker: about 0.16 m is gained relative to the case with low density snow.

Fig. 4. Ice thickness for some sensitivity experiments concerning snow cover. The solid black curve is for the case without any snow on the ice. The dashed curve shows the effect of keeping the snow density at 200 kg m⁻³. The reference curve corresponds to the simulation of Figure 2 (with a snow density of $500\;{\rm kg}\;{\rm m}^{ \hbox{-} 3}$). The red curve shows simulated ice thickness for a climate change experiment according to the RCP-2.8 Scenario (NCCS, 2018).

The relation between lake ice cover and climate change is an important issue. In several studies the effect of climate change on ice cover over larger lakes in mid- to high latitudes has been investigated (e.g. Brown and Duguay, Reference Brown and Duguay2011; Dibike and others, Reference Dibike, Prowse, Saloranta and Ahmed2011). These studies have made clear that lake ice cover is quite sensitive to warming, because freezing rates will be lower and the date of freeze-over shifts foreward. The assessment of the effect of warming on a small mountain lake is particularly complicated, because the evolution of the water depends on so many factors (e.g. Liston and Hall, Reference Liston and Hall1995). An early snowfall event in autumn, for instance, with melting of the falling snow in the lake, leads to lowering of the water temperature and facilitates early freezing. Snowfall shortly after the formation of ice will reduce the thickening rate of the ice cover. Changes of the inflow of water from the adjacent mountain terrain also have a strong impact on the cooling rate of the lake in autumn. In view of the the complexity of the thermal system of a mountain lake, any estimate of the effect of climate change on the lake ice cover will have large uncertainties. Nevertheless we have carried out a simple calculation for the projected wintertime temperature change and precipitation change in the year 2060 under the RCP-2.8 scenario (being + 1.6 K and 20% for the upper Engadin region, respectively; see NCCS, 2018). From a simple lake water temperature model (unpublished) we found that such a temperature increase would shift the date of freezing-over by about 11 days. This shift is in broad agreement with the values found by Dibike and others (Reference Dibike, Prowse, Saloranta and Ahmed2011) for comparable temperature changes. The red curve in Figure 4 shows that the impact is large: the ice thickness at day 90 would be reduced by about 30%. For other RCP scenarios with a larger temperature increase the changes would be even more dramatic. We believe that simulations over a longer time span are needed to analyse the climate sensitivity of lake ice cover in a statistically sound way. However, our simple calculation points to a large sensitivity of the ice cover of Lake St Moritz.

4. Summary and discussion

In this letter we have presented and tested a simple model for lake ice growth (SLIM) under alpine conditions. In contrast to more sophisticated models, we have chosen to express the upper temperature of the ice layer directly in terms of air temperature and snow cover. With freezing-over date, daily mean air temperatures, and snow depth on the ice as input, the model simulates the ice growth on Lake St. Moritz quite well (σ = 2.1 cm). This result was achieved by finding an optimal value for the thermal resistance of the snowpack on the ice. Applying the model to the nearby Lake Silvaplana also gave good results (σ = 3.8 cm), especially when the thermal resistance of the snowpack is increased (σ = 2.3 cm).

We emphasize that the conditions in the high alpine valley of the Engadin are relatively simple. In spite of the recent warming trends, rain events in winter are still very rare. The same applies to gales and blizzards.

We have been able to quantify the effect of grooming of the snowpack on ice growth. Apparently grooming is an efficient way to enhance the ice thickness. However, this point deserves further study. In the future we plan to carry out ice thickness measurements on groomed and non-groomed surfaces in a systematic way. Snow density measurements will also help to validate the approach we have taken.

Within the context of our simple approach, some model refinements can be made. The timescale τ, for instance, could be made a function of the ice thickness or snow depth. The meteorological input data could also be refined. However, given the small rms error of the simulated ice thickness, there appears not much to be gained.

We did not consider the situation in which the snow burden depresses the ice so much that water flows into the snowpack and snow-ice is formed. This may happen when the snow depth exceeds the critical value $h_{\rm s}^\ast$, given by

(8)$$h_{\rm s}^\ast{\approx} \displaystyle{{110} \over {{\bar{\rho }}_{\rm s}}}h_{\rm i}, \;$$

where $\bar{\rho }_{\rm s}$ is the mean density of the snow column in kg m⁻³. Once a water layer has been detected, the freezing rate of this layer can be obtained from the heat flow through the snowpack, which is estimated to be k _s(T _a − T _m). At the same time, the growth rate at the bottom of the original ice layer should be set to zero.

It would of course be very interesting to apply the model to other lakes in a similar setting, but the application of the model certainly has its limits. Large amounts of snow or frequent windy conditions will be problematical to handle. However, application of SLIM to other lakes in a similar setting could help to further assess the model and identify other possible improvements.