Brief overview of the Crocus snowpack model

Crocus is a one-dimensional, multilayer, physically based snowpack model that simulates in detail the time evolution of the physical properties of layers making up the snowpack, solely driven by meteorological conditions at the surface (air temperature, wind speed, relative humidity, incoming long-wave and shortwave radiation, snowfall and precipitation rate). For each numerical snow layer, i, within the snow-pack, the model prognostic variables are temperature, T(i) (K), dry density, ρ(i) (kg m⁻³), i.e. the density of the ice matrix, thickness, dz(i) (m), liquid water content, w(i) (kg m⁻³), snow age (days since snowfall) and grain-type properties. The latter are described by four empirical variables (dendricity, sphericity, size and metamorphic state) representing the snow microstructure. More details on the variables used to describe the snow layers in Crocus are given by Reference Bouilloud and MartinBouilloud and Martin (2006), Reference WillemetWillemet (2010) and Reference VionnetVionnet and others (2012). The energy budget of each snow layer, i, is written as follows:

where C_p (i) is the specific heat capacity of ice, depending on temperature, and L _f is the latent heat for ice fusion (3.33 × 10⁵ J kg⁻¹). On the right-hand side of the equations, W _p and W represent the liquid water input through rainfall or percolation, respectively. S _abs (i) represents the fraction of incoming shortwave radiation absorbed by layer i, L _net is the net longwave radiation flux, H and LE are the turbulent fluxes for sensible and latent heat, respectively, P is the sensible heat flux due to precipitation, Q _c (i) represents the divergence of the conduction flux within layer i and Q _g is the basal heat flux. Figure 1 provides an overview of the different processes accounted for by the Crocus snowpack model.

Fig. 1.

Overview of the processes accounted for by the Crocus-DEB snowpack model, showing the possible presence of a debris layer sandwiched between snow or ice below and a transient snowpack on top (adapted from Reference VionnetVionnet and others, 2012).

The computation of energy fluxes within and at the boundaries of the snowpack by Crocus has been described extensively in previous publications (e.g. Reference Brun, Martin, Simon, Gendre and ColéouBrun and others, 1989, Reference Brun, David, Sudul and Brunot1992; Reference Martin and LejeuneMartin and Lejeune, 1998; Reference WillemetWillemet, 2010; Reference VionnetVionnet and others, 2012) and is not repeated here.

From Crocus to Crocus-DEB

Debris-covered glaciers are modeled as a vertical stack of ice, debris and a potentially developing snowpack at the surface. In this regard, the basic principles of Crocus, i.e. the time evolution of the physical processes of a snowpack consisting of several numerical layers, remain unchanged. We describe below the handling of debris layers by Crocus-DEB, including the simplifying hypotheses made.

Variables. To differentiate between numerical layers consisting of debris and those of snow and ice, we use an arbitrarily out-of-range value for one of the prognostic variables used by Crocus to handle snow metamorphism (namely, the variable representing either sphericity or grain size). The only true prognostic variable considered for a debris layer is its temperature. In contrast to a snow or ice layer, the temperature of a debris layer can take values >273.15 K. Debris layers are attributed a given thickness, density, thermal conductivity and specific heat capacity. These properties remain constant throughout a given simulation. In this work, all the numerical layers making up the debris cover were attributed the same physical properties, but this is not imposed by the model structure. We further hypothesize that the liquid water content of debris layers is always zero, implying that liquid water percolation through the debris is instantaneous and that phase changes do not occur in the interstitial volume in the debris layers.

Processes. When there is a debris layer at the surface, the computation of the surface energy budget is carried out in a different way to the case where snow is present at the surface. When debris is present at the surface, different values than those for snow are used for the thermal emissivity, ε, roughness length, z ₀, and spectral albedo. The spectral absorption coefficient of the debris is considered infinite, i.e. all of the net shortwave radiation budget is affected to the uppermost debris layer. The thermal emissivity of the debris, ε _d, is kept constant. The spectral albedo of the debris is fixed throughout a given simulation. The sensible heat flux, H, is computed similarly to snow; however the debris can be attributed a different roughness length, z ₀. The latent heat flux, LE, is zero because the debris is considered to be always dry. Liquid water from rainfall is supplied to the uppermost debris layer. The corresponding heat flux is computed similarly to the rainon-snow case (i.e. rainfall temperature is taken equal to air temperature) (Reference Brun, Martin, Simon, Gendre and ColéouBrun and others, 1989, Reference Brun, David, Sudul and Brunot1992; Reference VionnetVionnet and others, 2012). For the case of snowfall on top of the debris (in a similar way to how Crocus handles snowfall on the ground) a snow cover may or may not develop, depending on the thermal state of the debris. When the debris temperature is below the melting point for ice, a snow cover develops on top of the debris. For higher debris temperatures, development of a snow cover depends on the balance between the heat content of the debris and the falling snow.

When debris is covered by snow, all surface fluxes are computed similarly to the standard Crocus simulation. We use the parameterization of snow albedo introduced by Reference LejeuneLejeune and others (2007), which serves to lower albedos, α _sd, of the snowpack under a given critical thickness:

where α _d (j) and α _s (j) are the albedo of debris and snow, respectively, in a given spectral band j. F_α is a weighting function defined as

where h _sd is the total snow depth over the debris, is the critical snowpack thickness under which the debris has an impact on the snowpack albedo () and x is an empirical coefficient (x = 0.33) (Reference LejeuneLejeune and others, 2007). This approach allows us to account for the fact that a thin snowpack on top of the debris cover or on top of the ground is influenced by the underlying absorbing surface, effectively lowering its albedo, and also represents the fact that thin snowpacks are often uneven and patchy, which tends to accelerate their melting.

Handling of numerical layers. Crocus features an evolved scheme to dynamically adapt the vertical grid mesh to the snowpack stratification (Reference Brun, David, Sudul and BrunotBrun and others, 1992; Reference WillemetWillemet, 2010; Reference VionnetVionnet and others, 2012), leading to automated splitting and aggregating of the numerical snow layers during a model run, depending on the similarity of their physical properties and on the vertical discretization of the snowpack. This feature is switched off for debris layers, i.e. the latter do not undergo any splitting/aggregation throughout a given model run, and keep a constant thickness. The thickness of the numerical debris layers has to be chosen sufficiently small to accurately solve the heat diffusion equation, with a minimum imposed by the numerical scheme and the time-step of the model. A good compromise was found to be a numerical debris-layer thickness of 1 cm and an internal model time-step of 15 min. In this work, the debris cover was split into as many 1 cm thick layers as necessary. Note that the thickness of each debris-cover layer can be set individually, and can be adapted depending on the application.

Below-debris temperature. The reference implementation of Crocus-DEB solves the heat diffusion equation (including the phase change in snow and ice) through the whole snow/ debris/ice continuum. However, for the purpose of comparing Crocus-DEB with previous models of debris-covered glaciers assuming that ice below the debris is at the melting point (e.g. Reference Nicholson and BennNicholson and Benn, 2006; Reference Reid and BrockReid and Brock, 2010), a modified version of Crocus-DEB was developed in which the temperature of sub-debris layers is forced to the melting point. At each time-step, the temperature of below-debris layers is raised to the melting point, and the corresponding amount of energy artificially added to the system is recorded. This model version is referred to as Crocus-DEB_m below, for the sake of brevity.

Field experiment

The research station Col de Porte (CDP) was used to test and evaluate Crocus-DEB. The CDP experimental snow and meteorological station has been recording the physical properties of the seasonal snowpack since 1959, along with the meteorological conditions. In addition to a detailed description of the site, Reference MorinMorin and others (2012) provide a quality-controlled snow and meteorological dataset spanning the period 1 August 1993 to 31 July 2011 at the site, at hourly time-steps, which includes the four components of the surface radiative balance, air temperature and relative humidity, wind speed and direction, precipitation (amount and phase), snow depth and snow water equivalent (SWE), and surface temperature. In addition, at weekly time resolution, in situ snow-pit observations are manually performed at three locations within the field site to measure the SWE, as well as the vertical profile of the physical properties of the snowpack (grain type, density, temperature, liquid water content). Due to differential shading and occasional snowdrift, the snowpack can be spatially heterogeneous (Reference MorinMorin and others, 2012). Nevertheless, the whole area is considered here as the reference debris-free test area (TA), TA₀, for our experiment.

Between February and April 2011, inside the CDP experimental field, a dedicated field experiment was carried out to investigate the impact of a well-characterized debris cover on the underlying snowpack, and to provide data to drive and evaluate the model. The seasonal snowpack was used as a surrogate for the ice mass underlying the debris in the case of a real debris-covered glacier. The gravels were deposited on top of the snowpack on 1 February 2011 at 12:00. One of the two test areas of 1.5 m × 3 m was covered with a 4 cm thick layer, the other with a 15 cm layer; these ‘debris-like’ layers consisted of rounded gravel (mixed rock types including igneous, sedimentary and metamorphic rocks extracted from the Isère river in a quarry near Grenoble, France) with a density of 1460 kg m⁻³ and a bead size ranging from 10 to 20 mm. These two test areas are referred to hereafter as TA₄ and TA₁₅, respectively, and their results are compared with those of the reference debris-free area, TA₀. Figure 2 provides an overview of the experimental set-up. Sensors were placed to continuously measure the temperature at the bottom, in the middle and at the top of artificial debris layers (Pt100 temperature probes), and the conductive heat flux at the bottom of the debris (Hukseflux HFP01). In addition, incoming and reflected shortwave radiation were measured on top of TA₁₅ with a pyranometer (Hukseflux NR01). From 8 February 2011 at 16:00 to 7 April 2011 at 08:00, approximately weekly visits were made to the study site to carry out manual measurements of snow depth and SWE above and below the artificial debris, as well as on TA₀. Manual measurements were carried out on one side of the TA, while automated measurements were carried out on the opposite side, to minimize disturbance to the continuous measurements. Manual snow corings providing the sum of SWE below and above the debris were performed on 8, 11 and 24 February and 21, 25 and 29 March 2011. Measurements limited to the upper snow layers above the debris were additionally carried out on 24 February and 11 and 15 March 2011.

Fig. 2.

Schematic view of the experimental set-up at Col de Porte. The two areas with artificial debris thicknesses of4 and 15 cm, referred to as TA₄ and TA₁₅, respectively, are displayed beside the debris-free area, denoted TA₀, which in practice corresponded to the regular snow observations area at CDP (Reference MorinMorin and others, 2012).

Numerical simulations

The Crocus-DEB model is driven using the same meteorological forcing as Crocus (Reference Brun, David, Sudul and BrunotBrun and others, 1992; Reference WillemetWillemet, 2010), i.e. hourly records of air temperature, T _a, relative humidity, H _a, incoming longwave and shortwave radiation, L _in and S _in, wind speed and snow- and rainfall. Model simulations were carried out on TA₀, TA₄ and TA₁₅, to numerically reproduce the conditions of the field experiment held at CDP. The values for roughness length, z _{0, s}, and thermal emissivity, ε_s, for snow were the usual values of 0.004 m and 1.0, respectively. The thermal emissivity, ε_d, thermal conductivity, k _d, and specific heat, C _p,d, for the debris were assigned the values 0.95, 0.7 W m⁻¹ K⁻¹ and 950 J kg⁻¹ K⁻¹, consistent with Reference Reid and BrockReid and Brock (2010) for their study on Miage glacier, Italy. The values of density, ρ _d, and albedo, α _d, of the debris were determined experimentally and the obtained values of 1460 kg m⁻³ and 0.2, respectively, were used for the numerical simulations. The value of the roughness length when debris is at the surface, z _0,d, was set to a value of 0.012, which is higher than that of snow but takes into account the fact that the test areas were surrounded by snow. Given that turbulent heat fluxes at the surface/atmosphere interface operate over spatial scales of a few metres at least, it was necessary to account for spatial heterogeneity of the surface in setting the z _0,d value. Nevertheless, results were not significantly modified when the z _0,d value was varied within a reasonable range. The main adjustable parameters used for the simulation at CDP are summarized in Table 1.

Table 1.

Values of the main adjustable parameters for the Crocus-DEB numerical simulation at Col de Porte and Changri Nup Glacier

The model runs were set up and initialized as follows. First, a full model run without any debris was carried out from 1 September 2010 at 00:00 to 12 April 2011 at 0:00, corresponding to the reference simulation on TA₀. In order to carry out the simulation including the presence of debris, on 1 February 2011 at 12:00, i.e. on the date when artificial debris was poured on top of the snowpack, the simulated vertical profile of the physical properties of the snowpack was altered by adding the required number of debris-cover layers to simulate a 4 and a 15 cm thick cover. The initial temperature of all debris layers was 273.15 K. Crocus-DEB model runs were then performed for TA₄ and TA₁₅ conditions starting on 1 February 2011 at 12:00, until 8 February 2011 at 16:00, the date of the first field measurement of the total snow depth and SWE below the artificial debris cover. The simulated profiles on 8 February at 16:00 were adjusted for TA₄ and TA₁₅, to match the observed snow depth and SWE below the artificial debris cover on each TA. Model runs for TA₄ and TA₁₅ were then initialized using the distinct starting conditions on 8 February at 16:00 and run until 12 April 2011. The same procedure was carried out when using Crocus-DEB_m.

Delineation of four typical periods

Based on the meteorological conditions and the presence or not of a transient snowpack on top of the debris, the measurement campaign was split into four periods (Figs 4–7; Table 2), and summarized along with their typical meteorological characteristics. Due to the peculiar conditions that prevailed during the 2010/11 winter, these periods were exceptionally homogeneous in terms of their meteorological conditions.

Table 2.

Surface state of debris-cover test areas and average meteorological conditions (air temperature, relative humidity, incoming longwave and shortwave radiation, and wind speed) during the periods P1, P2-1, P2-2 and P3 at CDP. Date format is dd/mm/yyyy

Snow on top of debris

Table 3 provides measured and simulated SWE and snow-depth data above the debris on three dates during the experimental campaign, after the two significant snowfalls on 20 February and 26–27 February (Figs 4c and 5c). It shows that these two integrated snow properties are well simulated by Crocus-DEB. These snowfalls are characterized by a SWE increase visible in Figures 6a and 7a.

Table 3.

Measured and simulated values for snow depth and SWE of the transient snowpacks overlying debris cover on TA₄ and TA₁₅ to

Total snow water equivalent

In terms of total SWE (i.e. taking into account snow below and above the debris layers), there is overall good agreement between the model results and the observations (Figs 6a and 7a). On TA₄, observed and calculated SWE agree well throughout the whole period of the experiment, except on 24 February, when the measured and simulated SWE are 145.5 and 114.7 kg m⁻², respectively. The agreement is very good in March, leading to a simulated complete melting of the snow cover on 30 March that is within 1 day of the total melt date identified using temperature measurements. On TA₁₅, in March, the melting rate is well reproduced by the model but there is systematic overestimation of the SWE (of 32, 19 and 23 kg m⁻² on 21, 25 and 29 March, respectively), leading to complete disappearance of the snow on 10 April, 6 days later than the field observations.

Albedo

Figure 7a shows the observed and simulated albedo on TA15. It shows overall agreement between the two, with lower values, ∼0.20, in the absence of snow, and higher values, up to 0.8, in the presence of snow at the surface. Consistent with other evaluations of Crocus, in the presence of snow the simulated albedo is generally higher than the observations (despite the implementation of the formula reducing the snow albedo for shallow snowpacks; see above). This may explain why the complete melting of the snow on top of the debris is delayed by a few days in the simulation (cf. the observations).

Heat flux and temperature within the debris

Heat fluxes and temperatures at the top, in the middle and at the bottom of the debris-cover layers are reasonably well simulated by the model for both debris thicknesses (Figs 6b–f and 7b–f). Indeed, the model satisfactorily captures the dynamics and range of variation of the temperature and heat flux within and at the boundary of the debris cover. The main discrepancies occurred on 25 and 26 February. For these days, on both TA₄ and TA₁₅, the model overestimates the melt rate of the thin snow layer occurring on top of the debris for just a few hours (total snowfall precipitation of only 10 kg m⁻²), explaining why temperatures and heat fluxes then show a large daily variability, in contrast to the observations. Consequently, the calculated melting during these two days is 20 and 10 kg m⁻² over TA₄ and TA₁₅, respectively, although it was insignificant from the measurements. Between 2 and 13 March, simulated temperatures were below the melting point, leading to slightly negative daily basal fluxes, although the debris always remained at or slightly below freezing. Consequently, the model erroneously simulated a cooling of the snow below the debris during period P2-2. Finally, the relative differences of the observed and calculated diurnal amplitudes of the debris temperatures and basal fluxes are higher over TA₁₅ than over TA₄ during period P1. This is also the case during period P3, but exclusively for debris temperatures and not for heat flux. Indeed, between 20 and 27 March anomalously low temperatures were simulated at the bottom of the debris over TA₁₅, albeit daily basal heat fluxes were well simulated. Although the date of the complete disappearance of the snow is well simulated over TA₄, the calculated cumulative heat flux on 29 March is 120 W d m⁻² higher than the observed cumulative flux, i.e. a discrepancy of ∼20%. In contrast, over TA₁₅, the simulated and measured cumulative fluxes are in good agreement, but the date of snow disappearance is simulated to be too late. These discrepancies clearly show that there are cases of error compensation in the model.

Snow covering the debris (period P2-2)

During period P2-2, i.e. when snow was covering debris on TA₄ and TA₁₅, field data confirm that the presence of snow had a particularly strong attenuating effect on heat transfer, as seen in Figure 8, where the snow layer was 10–30 cm thick. At TA₄ (Fig. 8a), temperatures at the three measurement levels (top, middle and bottom) always remain slightly below the melting point.T _{bottom–4–meas} ranged between 272.9 and 273.0 K with a daily mean of 273.0 K. Daily means of T _{mid–4–meas} and T _{top–4–meas} were 273.1 and 273.0 K, respectively, and their respective daily amplitudes of 0.1 and 0.2 K were barely higher than that of T _{bottom–4–meas}. At TA₁₅, debris always remained at the melting temperature (Fig. 8b) and the basal flux, Q _{c–bottom–15–meas}, was 0 (Fig. 8f). At TA₄ the daily basal flux was only −0.2 W m⁻² (Fig. 8e).

The results of the simulation during period P2-2 are given in Figure 8c, d, g and h. They show that the numerical simulations overestimate the diurnal amplitude of debris temperature. At both TA₄ and TA₁₅, simulated temperatures are always below the melting point, and minimal at midday. Temperature increases from the top of the debris layer toward the bottom. Although measured basal fluxes were always 0, or very close to 0, simulated fluxes are always negative and range between −1 and −3.5 W m⁻² at TA₄ and between 0 and −2.5 W m⁻² at TA₁₅. The model tends to simulate the downward propagation of a ‘cold wave’, although measurements indicate that the snow below the debris is almost perfectly insulated from the surface.

Melting period without snow on top of the debris (period P3)

During period P3, corresponding to the main melt period without snow on top of the debris, experimental results (Fig. 9a and b) indicate that the surface temperatures measured on the debris over both test areas (T _{top–4–meas} and T _{top–15–meas)} minimize and maximize at the same time of the day, at around 5.00 and 13.00 UTC. The diurnal amplitude of T _{top–4–meas} is 22 K, slightly lower than that of T _{top–15–meas} (26 K). At the surface, the debris follows a thermal cycle, which appears to be similar regardless of the debris thickness, even if the maximum temperature is 5 K higher over TA₁₅ than over TA₄. In the middle of the debris layer, the diurnal amplitude of the temperature measured at -2 cm on TA₄ (T _{mid–4–meas}) is 15 K, three times higher than at −7.5 cm on TA₁₅ (T _{mid–15–meas}) (Fig. 9a and b). Also visible is a clear delay for T _{mid–15–meas} compared to T _{mid–4–meas}, the minimum and maximum of T _{mid–15–meas} being reached at 9.00 and 16.00 UTC, respectively, ∼4 and 3 hours after T _{mid–4–meas}, respectively. The heat flux propagating within the debris is then all the more attenuated and delayed as the debris is deeper in the layer. At the bottom of the debris layer, sensors were installed at the snow/debris interface, preventing the temperature from exceeding the melting point. At the bottom of TA₁₅, the temperature, T _{bottom–15–meas}, always remains at 273.15 K (melting), although on TA₄, underlying snow melts for only 6 hours a day, between 10.00 and 16.00 UTC. The rest of the time, T _{bottom–4–meas} is below the melting point, reaching its minimum value at 6.00 UTC (271.3 K). The underlying snow thus alternates between thawing and refreezing along the diurnal cycle.

The measured conductive heat fluxes at the snow/debris interface on TA₄ and TA₁₅, Q _{c–bottom–4–meas} and Q _{c–bottom–15–meas}, respectively, are plotted in Figure 9e and f. On TA₄, Q _{c–bottom–4–meas} is positive during daytime between 8:00 and 18:00, with a maximum value at 14:00 (100 W m⁻²), and remains slightly negative (∼−7 W m⁻²) the rest of the time. On TA₁₅, Q _{c–bottom–15–meas} always remains positive, with significant values only between 10:00 and 24:00, a maximum of 50 W m⁻² being reached at around 16:00. Even if Q _{c–bottom–4–meas} is negative at night, which means that refreezing occurs under TA₄ debris, its daily mean is 19 W m⁻² (cf. 14 W m⁻² for Q _{c–bottom–15–meas}), indicating that daily melting is higher below TA₄ than below TA₁₅.

Figure 9c, d, g and h show the results of the simulations. Comparison of Figure 9a and c forTA₄, or Figure 9b and d for TA₁₅, shows that simulated temperatures agree fairly well with observations, as also confirmed by correlation coefficients in Table 4 (r ² values 0.74–0.95, except for temperature at the bottom of TA₁₅ where r ² values are only 0.41 but the temperature variations are small). Nevertheless, the model fails to consistently reproduce the maximal heating at the surface of TA₄ and TA₁₅, with underestimations as high as 7 and 9 K, respectively. At the surface and in the middle of TA₄, mean daily values of T _{top–4–sim} and T _{mid–4–sim} are 1.3 and 0.6 K, respectively, lower than the measurements (Table 4). At TA₁₅, daily T _{top–15–sim} is 1.3 K lower than T _{top–15–meas} but daily T _{mid–15–sim} is 0.3 K higher than T _{mid–15–meas}. The diurnal amplitude of T _{mid–4–sim} is also weaker than that of T _{mid–4–meas} (Fig. 9a and c), although that of T _{mid–15–sim} is higher than that of T _{mid–15–meas} (Fig 9b and d). Finally, at the bottom of both debris layers, daily mean simulated temperatures, T _{bottom–4–sim} and T _{bottom–15–sim}, are slightly lower (∼0.2 K) than observations. There is good agreement between T _{bottom–4–sim} and T _{bottom–4–meas} with observed and simulated melting conditions encountered between 10.00 and 16.00 UTC, but T _{bottom–15–sim} remains slightly below the melting point between 4.00 and 12.00 UTC, although observations show that underlying snow is still melting.

The comparisons between Figure 9e and g and between Figure 9f and h show that the model is able to consistently reproduce the timing of the daily cycle of the heat flux at the bottom of TA₄ or TA₁₅ (with r ² between simulated and measured fluxes as high as 0.97 and 0.88, respectively; Table 4), but fails to properly simulate its daily amplitude. Indeed, the daily amplitude of Q _{c–bottom–4–sim} (220 W m⁻²) is twice as high as that of Q _{c–bottom–4–meas} and the daily amplitude of Q _{c–bottom–15–sim} (75 W m⁻²) is one-third higher than that of Q _{c–bottom–15–meas.} Looking at daily means, the agreement is much better, with daily values of Q _{c–bottom–4–sim} and Q _{c–bottom–4–meas} equal to 17 and 19 W m⁻², respectively, and daily values of Q _{c–bottom–15–sim} and Q _{c–bottom–15–meas} equal to 3 and 14 W m⁻², respectively.