An Energy and Mass Model of Snow Cover Suitable for Operational Avalanche Forecasting

Abstract

A numerical model has been developed to simulate energy and mass evolution of snow cover at a given location, as a function of meteorological conditions: precipitation, air temperature, humidity, wind velocity, and incoming short-wave and long-wave radiation.

This model, named CROCUS, was first tested on a well-instrumented field during a whole winter, showing its ability to simulate the important phenomena affecting the evolution of the snow layers: high temperature gradients, wetting, compaction, and melting-freezing cycles. A second test was conducted at two locations in the French network used for operational avalanche forecasting. Though the weather observations are made there only twice daily, the snow profiles calculated by the model were very close to those obtained once a week by a pit observation. CROCUS proved itself sufficient to be considered now as a useful objective tool for operational avalanche forecasting.

Introduction

The mechanical equilibrium of snow cover, at any given moment, is governed by the thermal and morphological state of its different layers: temperature or liquid water content, density, snow-grain type and size. This state, at a given location, is governed by the meteorological conditions which prevailed since the beginnning of the snow coverage.

In France, the avalanche-forecasting services collect information once a week from pits on the internal snow-cover structure and the forecaster assumes the evolution of the snow cover between two pit dates from the daily meteorological observations. But this step is difficult and subjective, since many parameters having considerable spatial variability work together: slope, orientation, air temperature, wind velocity, humidity, cloudiness, and snow albedo. The forecaster needs the help of an objective tool.

Knowledge of snow properties now allows the development of a numerical simulation of energy and mass evolution of the snow cover as a function of past and present weather conditions. Such a model would supply relevant information for the day-to-day avalanche forecasting. In past years, various evolution models have been developed, but principally for hydrological applications (Obled, 1971; Anderson, 1976; Morris and Godfrey, 1978). Navarre (1975) developed the first model suitable for avalanche forecasting.

The purpose of this paper is to describe a new energy and mass model of the snow cover which is well suited to monitoring the evolution of its different natural layers. This model was first tested at the measurement site at Col de Porte (French Alps) where complete meteorological observations are continuously and automatically recorded. As a second step, it was tested at two locations selected from the French snow-weather network where the meteorological conditions were measured twice daily. During the two experiments, the results of the model were compared with the observations made within the snow-pack. It allowed an estimation of the quality of the model which appeared sufficiently high to consider this model as a useful and necessary objective tool for day-to-day avalanche forecasting.

Energy and Mass Balance of the Snow Cover

The snow-pack evolution models developed for hydrological applications were principally developed to describe snow-cover melting. Anderson (1976) took into account all the phenomena affecting the energy and mass balance, and his model was efficient for the simulation of bottom-water run-off. Morris and Godfrey (1978) described the internal water transport by a physical equation but did not introduce the internal absorption of solar radiation which is then supplied to the snow surface. When the snow surface melts, it does not change the energy balance but, in the opposite case, it warms the snow surface more than in reality and the energy balance may be affected. In some cases, such a simplification cannot describe snow melting when it occurs a few centimeters below the surface which is still at a negative temperature.

These previous models are efficient for hydrological applications but the numerical layers they use do not match the natural layers of the snow-pack and they are insufficiently suited to following the energy history of each natural layer, which is the main information for characterizing its mechanical properties.

This heading describes the choice we have made to take into account all the phenomena affecting the evolution of the snow cover. To simulate a snow-pack, we consider it as unidimensional. The energy exchanges are projected perpendicular to the slope and are positive when they are supplied towards the snow-pack.

Long-wave radiation: Ql

The radiative properties of snow in the range 5–40 µn imply radiative exchanges strictly confined to the snow surface.

(1)

where Q↓ is the incident atmospheric long-wave radiation, T_s is snow-surface temperature, σ is the Stephan constant, and ε_s is snow emissivity. ε_s depends on the wavelength and on the angle of incidence. It generally varies from 0.98 to 1. We chose it equal to 1.

Short-wave radiation: Qs

Incoming solar radiation is partly reflected by snow. The remaining part penetrates through the snow where it is gradually absorbed. Snow reflectance depends strongly on wavelength, grain-size, and impurity content. As the solar spectrum distribution varies with cloudiness, snow albedo α is not constant for a given snow layer. Absorption β depends on wavelength, grain-size, density, and impurities (Bohren and Barkstrom, 1972; Sergent and others, 1987). Since the penetration of solar radiation with depth is an exponential function e ^−βz with β varying strongly with wavelength, the solar_ absorption cannot be described by β unique mean value β. Therefore, the solar spectrum was divided into three spectral bands, 0.3–0.8, 0.8–1.5, and 1.5–2.8 μm, over which α and β are considered constant and depend only on the grain-size and density according to Sergent (unpublished) and Bohren and Barkstrom (1972).

At a depth z below the snow surface, the solar flux Qs is given by:

where Rs_i is the incoming solar radiation in the spectral band i. The radiation absorbed by a layer of thickness ▵z at a depth z is equal to:

Turbulent exchanges between the snow surface and the atmosphere: Qh and Qe

Deardorff (1968) expressed the sensible turbulent flux Qh and the latent turbulent flux Qe as:

where p_a is air density, Cp_a is specific heat of air, U_a is wind velocity at a given height above the snow surface, T_a is air temperature at the same level, Ls is the ice latent heat of sublimation, P_a is the atmospheric pressure, M_vM_a is the ratio between water-vapor and dry-air molecular weights, E_i(T) is the saturation vapor pressure above a flat ice surface at the temperature T, and T_s is snow-surface temperature. Ch and Ce are turbulent transfer coefficients which are generally assumed equal. They depend on snow-roughness length, on the height where T_a and U_a are measured, and on atmospheric surface boundary-layer stability. Above a snow-covered surface, the air temperature is generally higher than the snow-surface temperature in such a way that the boundary layer is stable. In this case Ch and Ce are very low (Deardoff, 1968) and, under slight wind conditions, the energy supplied to the snow by heat conduction through the air and by vapor diffusion due to vapor gradients in the air may be higher than by turbulent transfer. So, we prefer to suggest the following expression for Qh and Qe:

where a and b must be adjusted experimentally for a given location. In the above equations Ls must be replaced by (Ls - Lw) and E_i(T) by E_w(T₀) when T_s is equal to 0°C, since possible condensation would be liquid instead of solid when T_s equal to the melting point. Lw is the water latent heat of melting and Ew(T₀) is vapor pressure over a flat water surface at the melting-point temperature T₀.

Heat exchanges due to precipitation: Qr

Snow is supposed to fall at the snow-surface temperature and rain at the air temperature. When it is raining, liquid water is introduced within the upper snow layer at the temperature T₀ and the energy Qr equal to MrCp_w(T_a - T₀) is supplied to the upper layer. Cp_w is the water specific heat at T₀, and Mr is the liquid water mass.

Heat conduction through the snow-pack: Qc

where λ is the effective snow-conduction coefficient. Yen (1981) showed that most experimental λ measurements may be described by the following formula:

where λ_i is the ice conduction coefficient, p_n and p_w are snow and liquid water density. Accurate measurements of λ during experiments on thermal convection in snow confirmed this formula (Brun and Touvier, 1987). This coefficient is called effective since it includes heat fluxes due to vapor diffusion through the snow-pack and which may be formally considered as thermal conduction:

where ρ_v is vapor density, and D is the vapor diffusion coefficient in snow. Therefore, Qc must be limited by Qv in Yen’s formula.

Water movement through the snow-pack

When a snow layer is wet, liquid water may flow downward. According to Colbeck (1972), water flow occurs when water saturation exceeds the irreducible water saturation. The vertical water flux U_w may be expressed as:

where g is the gravity acceleration, U_w is water viscosity at 0°C, k is intrinsic snow permeability, Sw is water saturation, and Swi is the irreducible water saturation. Swi depends on snow type and grain-size, and is typically around 9–10% per mass (Denoth and others, 1979; Brun, 1989). Water run-off at the snow-pack bottom is supposed to penetrate integrally through the ground.

Snow settling

The snow layers settle by the combined effect of grain metamorphism and the weight of the upper layers. The latter effect may be described by the mean viscosity. A settling law was established by Navarre (1975):

where e is the layer thickness, σ is the vertical stress, di is the time interval, and f(d) is a function of the snow type. In the present case, our model does not describe snow-grain evolution and we therefore considered f(d) constant and equal to 0.4. Such a simplification may have serious consequences regarding the settling of depth-hoar layers, but no satisfactory solution can be found without explicit modeling of snow-grain metamorphism.

Heat transfer between snow and ground: Qg

Because of the inter-annual variations of ground temperature, an energy Qg, generally positive, is supplied to the bottom of the snow-pack. It depends on climatological and pedological conditions encountered at the specific location. Generally, Qg decreases slowly during the winter accumulation period and decreases significantly during the melting period when cold water flows through the ground. However, we consider Qg constant in the model, depending only on the location.

The above considerations about snow energy and mass balance show that the following input data must be available to simulate the snow-pack evolution: liquid and solid precipitation, wind velocity, air temperature and humidity, short-wave and long-wave incoming radiation, and thermal flux supplied from the ground.

Numerical Solutions

To calculate its temperature, density, and liquid-water content profiles, the snow-pack is divided into layers parallel to the slope. Energy transfers are projected perpendicular to the slope. The thickness of a layer ;’, dz(i), is variable versus depth and time. Since the variations of greatest amplitude occur at the snow surface, the thickness of the upper layers is smaller than the thickness of the bottom layers.

Heat conservation in an internal layer i may be written as:

where Cw(i) is the liquid-water mass within the layer considered, W is the liquid-water balance of the layer due to percolation, and Cp is the ice specific heat which is a linear function of the temperature T(i).

For the surface layer i, the equation becomes:

where Qc is the heat flux due to conduction with the layer just below, and W is the water due to precipitation.

For the bottom layer, the conservation equation becomes:

Because of possible phase changes between ice and liquid water, the energy-conservation equations cannot be solved simultaneously to calculate both variables T(i) and Cw(i) at time (t + dt) from their values at time t. We have chosen to calculate first the temperature T(i) of each layer i at time (t + dl) involved by the energy equation of the layer and then to impose phase changes necessary to have liquid water only when the temperature is equal to the melting point T₀. T₀ is assumed to remain constant and equal to 273.16°K.

To calculate the temperature of each layer, we use the classical resolution technique of Cranck and Nicholson which is an implicit method using finite differences. It requires the linearization of each heat exchange with regard to the temperature. The implicit method is necessary to get good accuracy for the exchanges near the surface. In our case, it is significant because Qh, Qe, and Ql depend strongly on the snow-surface temperature T_s. When T_s is equal to the melting point T₀ at time t, a preliminary computation of the surface-energy balance is made to determine whetherT_s will remain equal to T₀ at time (t + dt). In this case, the computation of the surface-energy balance is exact during the whole time interval and the implicit method is not used during that interval, since it would lead to a value of T_s strictly greater than T₀, so affecting strongly the heat exchanges.

After calculation of the temperature profile at time (t + dt), the model makes the necessary phase changes corresponding to the possible freezing of wet layers or their possible melting. The computation of water percolation is then made after introducing possible rain water. It uses a linearization of the percolation equation to avoid any numerical instability.

Settling is then taken into account by decreasing the thickness of each layer corresponding to the increase of density. Possible new snow layers are added to the snow-pack.

Before using the new temperature, density, and liquid-water content profiles for a new run, each layer depth is tested to avoid layer thicknesses below 0.5 cm and to avoid a number of layers greater than 50. In such a case, adjacent layers are combined according to the following rules, if possible: no combination of layers derived from two different snowfalls, combination of adjacent layers of closest density, keeping the thickness of the 15 upper layers below 1 cm. Such rules allow matching of the numerical layers to the natural layers of the snow-pack. This allows us to obtain, for a given layer, its complete “energy history” since its formation and so to predict its morphological and mechanical characteristics with regard to the snow-pack stability (high temperature gradients, wetness, compaction, and melting-freezing cycles). Because of the possibility of numerical layer thickness reaching 0.5 cm, the time interval was chosen equal to 15 min.

Test of the Model on a Well-Instrumented Site

To test the model, measured snow profiles over a whole snow season must be compared with the profiles calculated by the model. Possible differences may also be involved due to defects in the model, such as a poor knowledge of the prevailing meteorological conditions which are used as input data during the simulation. Therefore, a well-instrumented location was chosen for the first test of the model.

We chose the measurement site of our laboratory at Le Col de Forte, located in the Massif de la Chartreuse in the northern French Alps, 1320 m.a.s.l. The site is located in a slightly windy glade. Mean annual precipitation is around 2 m. Typical climatic conditions may be characterized by sequences of cold and warm periods during which rainfall may be seen even during the coldest months. Continuous snow cover usually lies from late November to the beginning of May. Deep snow layers are wet most of the time and the upper snow layers are submitted to varied conditions depending on the weather. Because of its variable climatic conditions, the chosen location is well suited to testing how the model is able to simulate the response to the main phenomena affecting the snow-pack: wetting, partial or total refreezing, rapid accumulation, and rain on snow.

Input data measurements

The following parameters necessary to the simulation were measured every hour: air temperature and humidity, liquid and solid precipitation, wind velocity, net solar radiation, and long-wave incoming radiation.

Snow-pack profile measurements

Hourly snow-temperature profiles were measured automatically by the following original device: after snowfall, a small plate was placed on the snow surface. This plate slides along an electric wire hanging from a portico. A small mercury switch connected to the plate makes an electrical contact with the wire. The plate carries a small platinum thermistor. Since no artificial vertical heat transfers are involved by the device, it allows us to obtain very accurate snow-temperature profiles (accuracy of ±0.02°C) and to measure accurately the settling of each layer (accuracy of ±2 mm). Snow depth, snow-surface temperature by long-wave output radiation, and bottom-water run-off were measured hourly.

Once a week, snow temperature, density, and liquid-water content profiles were determined on the site by a classical snow pit.

Simulation results

Continuous snow cover remained from 17 December 1986 until 27 April 1987. Figure 1 describes the observed snow depth and layering. Figure 2 describes the observed air temperature, precipitation, and bottom run-off.

Partial tests

Each of the physical laws introduced into the model was first tested separately choosing for each of them suitable periods.

The thermal conduction scheme was tested by simulation of the internal temperature at the given level of a plate using as boundary conditions the measured temperatures at the level of the plate just below and just above. The levels were chosen in order to have only heat conduction as energy transfer, i.e. the snow layers were dry and far enough below the surface to neglect the incoming solar radiation. Figure 3 describes the results of this simulation during 8 d.

Fig.1. Snow depth and compaction of each layer observed during the 1986–87 winter at Col de Porte. Missing snow-depth data are shown by the dashed line.

Fig.2. Maximum and minimum daily temperature, snow and rain precipitation, and bottom-water run-off observed during the 1986–87 winter at Col de Porte.

Fig.3. Comparison between measured (–) and simulated (–-) internal temperature during an S d run to test the conduction scheme.

The coefficients a and b used to calculate Qh and Qe were fitted from snow-surface temperature measurements during a cloudless and windy night in February.

Cold and sunny days in January were used to test the absorption of incoming solar radiation. The melting period was used to test the percolation scheme.

Global test

After having checked each parameterization, a complete simulation was conducted during three periods.

The first period lasted from 17 December 1986 until 25 January 1987. It was interrupted because of a breakdown of the data logger. The initial profile was determined when the snow depth was only 22 cm. Figure 4 shows a comparison between the measured and the observed snow depth. A small under-evaluation increasing versus time may be seen. Since the same defect was observed simultaneously concerning the snow-pack water equivalent, even during cold periods when no melting was either observed or simulated, we concluded it came from an under-evaluation of the snow precipitation, due to wind or evaporation effects in the rain gauge. Despite this defect, we noted a good agreement, principally for describing the settling of fresh snow.

Fig.4. Measured (––) and simulaled ( —) snow depth during the first test period al Col de Porte.

In Figure 5 are plotted the measured and simulated snow-surface temperatures during a cloudy and cold period followed by a warm and sunny one. Despite the variation of meteorological conditions encountered during this test period, the correlation is equal to 0.95 and the absolute error is equal to 2.1°C. The very good agreement validates the parameterizations used in the model for describing the energy exchanges between the snow surface and the atmosphere.

Fig.5. Measured (––) and simulated (—) snow-surface temperature during a 10 d run from the first period.

Figure 6 compares measured and simulated snow profiles after 34 d of simulation. Temperature profiles are very similar and snow depths show only a 10 cm difference.

Fig.6. Comparison between the simulated snow-temperature profile and the measured temperature (o) made at the level of the plates after a 34 d simulation.

The second period lasted 56 d from 4 February until 3 March 1987. It was interrupted when simulation results and observations were considered to be too different.

Figure 7 shows the measured and simulated snow depth during this period. As during the first test period, we remark that the difference increases versus time and is explained by an underestimation of the snow precipitation. However, the similar shape of the two curves between snowfall occurrences expresses a good description of energy exchanges and settling.

Fig.7. Measured (––) and simulated (—) snow depth during the second test period at Col de Porte.

The third period from 9 April until 26 April 1987 was the melting period. Only two slight snowfalls were observed and the model ization shows perfect agreement between the observed and the measured snow depth as seen in Figure 8. Figure 9 compares the observed and simulated bottom-water run-off. A systematic difference may be seen but the shapes of the two curves are very similar. The difference is due to weaker incoming solar radiation around the lysimeter than around the radiation sensors because of the proximity of a forest.

Fig.8. Measured (––) and simulated (—) snow depth during the third test period at Col de Porte.

Fig.9. Measured (––) and simulated (—) water run-off during the third test period at Col de Porte.

The simulation results during the whole season show that the model is efficient for simulating the different phenomena affecting the snow cover and its stability: surface snow-layer wetting, superficial or deep refreezing, layer settlement, and water run-off. All these events are governed by the snow-cover energy balance which depends strongly on the snow-surface temperature, in such a way that only a complete modelization of heat exchanges at the surface and within the snow cover allows determination of their occurrence.

Figure 10 shows the daily value of each type of heat exchange during the winter. Great variability may be noted, showing the interest of the model as an objective tool for the avalanche forecaster to quantify how the weather conditions affect the snow-pack.

Fig.10. Energy contribution of the different heat exchanges involved in the snow-pack energy balance depending on weather conditions at Col de Porte during the three test periods.

Test of the Model Within the Operational French Avalanche Forecasting Framework

The good quality of the model at Col de Porte is partly due to the good quality of the meteorological data used as input data during the simulation. Such data are not so easily available within the framework of operational avalanche forecasting, which is concerned with a large area. In the French Alps, the operational network has around 80 snow-weather observation points which produce two daily SYNOP-like weather observations, at 08.00 and 13.00 h. Air temperature, wind velocity, precipitation, and cloudiness are measured there.

To test the model, two locations were chosen where different climatic conditions usually prevail. The first one is at PAlpe d’Huez, Isère, 1900 m a.s.l. We chose a 9° south-facing slope. The second one was at Le Monestier, Hautes Alpes, 2100 m a.s.l., where we chose a 2° north-facing slope with high mountains southward, hiding the sun for an important part of the day. Empirical functions were established to extrapolate hourly values of T_a , V, and Q

from the two daily observations. Air humidity was deduced from measurements made at nearby locations. Incoming short-wave radiation was calculated from topography and an empirical function was used to evaluate the absorption by the clouds. Snow albedo was estimated from the observed snow-grain type.

The model results after simulation over a week were compared with the profiles of temperature, density, and sometimes liquid-water content measured once a week by an observer from pit. This was then used as the initial profile for the next weekly simulation. Such an approach fits the actual behaviour of the avalanche forecaster, who gets snow profiles once a week from each snow-weather station and guesses how these profiles are supposed to evolve with regard to the meteorological conditions.

The results obtained over the whole 1987–88 season at Le Monestier are shown in Figure, 11. Observed and simulated profiles are always in good agreement and the differences are of the same order as the natural spatial variation. Sometimes, near the surface, temperature profiles are different when, during cold and sunny days, the model computes in fresh surface snow layers a temperature profile characterized by a relative maximum a few centimeters below the surface, while the surface remains cold due to the long-wave radiation. Such a profile may be seen in Figure 11 for 26 January. They were sometimes observed but with such local gradients it is difficult to measure the temperature profile accurately. However, the main characteristics of the snow-pack are well simulated over the whole season. It remains dry until the end of March and the calculated temperature gradient is very close to the observed one. This information is very important for dry snow-packs, since temperature gradient governs dry-snow metamorphism, depending on its value in respect to thresholds to allow rounded crystals, faceted crystals, or depth hoar to grow. During late March, melting-freezing cycles occurred at the surface but liquid water did not flow through the whole snow-pack until the middle of April. Simulation and observation are again in very good agreement during this wetting period.

The results were of similar quality at the Alpe d’Huez observation point. However, despite a good global behaviour of the model, some defects were noted during the simulation. Since the viscosity used in the model does not take into account the snow type, the model calculates too great a density in the depth-hoar layers at the bottom of the Le Monsestier snow-pack. There are few values available in the literature to describe this viscosity. The second defect concerns water percolation through the snow layers, which is assumed to be unidimensional in the model. In Nature, water flow occurs preferentially within channels (Colbeck, 1975) and it cannot be easily taken into account in a mathematical one-dimensional model. This explains the difference found between the last observed and calculated liquid-water content profile at Le Monestier as seen in Figure 11. The two profiles are different, but the model and the observation show a snow-pack which is also wet at the surface as at depth, which is the important information in such a case.

Fig.11. Comparison between the snow-temperature and liquid-water content profiles simulated and observed from a pit once a week during the whole season at Le Monestier.

Conclusion

The two test experiments have shown that the model named CROCUS was efficient for simulating the snow profiles from the prevailing meteorological conditions. When only poor meteorological data are available, the loss in quality involved is not too important, in such a way that the model simulates well the different phenomena affecting the snow-pack: temperature gradient, surface wetting, refreezing of the wet-snow layers, melting, and bottom-water run-off. So, it may already be considered as a valuable tool for operational avalanche forecasting. However, some problems remain, concerning principally the settling of particular snow layers. This means that to improve this model simulation of snow metamorphism to determine the snow type of each layer is necessary.

Acknowledgements

The authors thank the people of the Centre d’Études de Ia Neige who were involved in the instrumentation of the Col de Porte laboratory, the avalanche forecasters of the Centre d’Études de la Neige who ran the model during the operational test, and also the ski resorts of l’Alpe d’Huez and Le Monestier which supplied the observations necessary for the test.