The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2

Abstract. Detailed studies of snow cover processes require models that offer a fine description of the snow cover properties. The detailed snowpack model Crocus is such a scheme, and has been run operationally for avalanche forecasting over the French mountains for 20 yr. It is also used for climate or hydrological studies. To extend its potential applications, Crocus has been recently integrated within the framework of the externalized surface module SURFEX. SURFEX computes the exchanges of energy and mass between different types of surface and the atmosphere. It includes in particular the land surface scheme ISBA (Interactions between Soil, Biosphere, and Atmosphere). It allows Crocus to be run either in stand-alone mode, using a time series of forcing meteorological data or in fully coupled mode (explicit or fully implicit numerics) with atmospheric models ranging from meso-scale models to general circulation models. This approach also ensures a full coupling between the snow cover and the soil beneath. Several applications of this new simulation platform are presented. They range from a 1-D stand-alone simulation (Col de Porte, France) to fully-distributed simulations in complex terrain over a whole mountain range (Massif des Grandes Rousses, France), or in coupled mode such as a surface energy balance and boundary layer simulation over the East Antarctic Ice Sheet (Dome C).


Introduction
Simulating the time and space evolution of the snowpack is key to many scientific and socio-economic applications, such as weather, hydrological (flood predictions and hydropower) and avalanche risk forecasting in snow-covered areas (Armstrong and Brun, 2008). When snow is present on the ground, it drives profound changes to all fluxes taking place at the interface between the Earth's surface and its atmosphere. Within the cryosphere, the seasonal snowpack is a very significant climate forcing , with a major impact on the energy budget of the soil and the atmosphere. At present, three major classes of snowpack models are used for various applications (Armstrong and Brun, 2008): single-layer snow scheme, scheme of intermediate complexity and detailed snowpack models. The main differences pertain to the description and the parameterization of the properties of the interior of the snowpack and the associated processes.
Snowpack models of the first class are generally included in numerical weather prediction (NWP) and climate models. In such models, the snowpack is represented as a single ephemeral soil layer featuring specific properties, such as a high albedo, a low thermal capacity and a low thermal conductivity. The snowpack is often represented with a fixed density. At present, despite major flaws in the quality of their representation of the physical properties of snow (Etchevers et al., 2004), they are commonly used in numerical weather prediction (NWP) and global climate models (GCM) (Douville et al., 1995) since they are relatively inexpensive, have relatively few parameters, and capture first order processes. Two snow schemes of this kind (D95: Douville et al., 1995, EBA: Bazile et al., 2002 are currently implemented in SURFEX (Le Moigne et al., 2009;Salgado and Le Moigne, 2010), within the Interactions between Soil, Biosphere, and Atmosphere (ISBA) land surface model (Noilhan and Planton, 1989), and are used in the operational NWP and  Acknowledging the limitations of single-layer schemes, snowpack schemes of intermediate complexity were developed to account for some internal processes such as snow settling, water percolation and refreezing. These schemes generally vertically discretize the snowpack with a prescribed number of layers (from 2 to 5, generally) (Boone and Etchevers, 2001;Loth and Graf, 1998;Lynch-Stieglitz, 1994). In these schemes, most snowpack physical properties are parameterized as a function of snow density, which is a surrogate for taking into account snow ageing (Boone and Etchevers, 2001). A snow scheme of this kind, named ISBA-Explicit Snow (ES), is currently implemented in SURFEX, within the ISBA land surface model (Noilhan and Planton, 1989;Boone and Etchevers, 2001), and is used operationally for hydrological applications in Météo-France (Habets et al., 2008). Many intermediate complexity snowpacks schemes exist, such as JULES (Best et al., 2011), CLASS (Brown et al., 2006), the Community Land-surface Model (CLM) (Oleson et al., 2010), WEB-DHM (Shrestha et al., 2010), and Snow 17 (Anderson, 1976). Models of this kind have been recently implemented within NWP and Earth's system models such as HTESSEL (Dutra et al., 2010) and RACMO (Kuipers Munneke et al., 2011).
Finally, a few detailed snowpack models belong to the third class and account explicitly for the layering of its physical properties. They include a more or less explicit description of the time evolution of the snow microstructure. This includes the models SNTHERM (Jordan, 1991), Crocus (Brun et al., , 1992 and SNOWPACK (Bartelt and Lehning, 2002). The representation of the grain morphology developed for Crocus and later implemented in SNOWPACK is based on semi-quantitative notions such as the dendricity and sphericity of snow grains, which can only be quantified using demanding image analysis processing (Lesaffre et al., 1998). Nevertheless, such models are best suited for reproducing the evolution of a snow season under the forcing of meteorological conditions, as demonstrated by the results of the Snow Model Intercomparison Project (Etchevers et al., 2004). Operationally, they are used in the field of avalanche risk forecasting, where the knowledge of detailed information on the vertical layering of the snowpack is critical (Durand et al., 1999;Rousselot et al., 2010). Regional or global simulations in coupled mode have been seldom carried out due to high computational costs .
Since its initial development, the snowpack model Crocus has been used in a stand-alone mode or coupled with various land surface models in a variety of environmental contexts. Some of the corresponding studies have constituted major scientific leaps in terms of the development and use of snowpack models. Indeed, Crocus has been the first model to simulate the metamorphism and layering of the snowpack (Brun et al., 1992). It made possible the first realtime distributed simulation of the snowpack over an alpine region for operational avalanche forecasting (Durand et al., 1999). In the 1990s, Crocus has been extensively used for the first physically-based studies to assess the impact of climate change on alpine snow climatology  and river discharges (Braun et al., 1994 of Crocus were implemented in the land-surface scheme of the regional climate model MAR to study snow/atmosphere interactions in polar regions (Gallée et al., 2001). Model limitations have also been highlighted. They concern mainly the interactions of the snowpack with its environment. In the first version of Crocus , the conductive heat flux at the snow/soil interface was set at a typical value observed at the experimental site of Col de Porte (1325 m altitude, French Alps). Several studies showed that this assumption fails under different climate or snow conditions: interaction between road surface and the overlying snowpack (Bouilloud and Martin, 2006), subarctic snowpack (Jacobi et al., 2010) or snowpack over a tropical glacier (Lejeune et al., 2007;Wagnon et al., 2009). To overcome this limitation, Crocus was coupled to ISBA by Bouilloud and Martin (2006) and this coupled version was further used to study the mass balance of the moraines over a tropical glacier (Lejeune et al., 2007), and for an intercomparison with several other snowpack models in terms of SWE (Snow Water Equivalent) simulations in Southern Quebec (Langlois et al., 2009). However the further development and use of this coupled version was not pursued and it is now obsolete. Crocus also did not include a representation of the snow-vegetation interaction which is crucial to simulate properly the snowpack evolution in forested areas (Rutter et al., 2009). Finally, only a limited number of studies refer to direct coupling of Crocus with an atmospheric model Durand et al., 2005). This article presents the current status of the snowpack scheme Crocus, now that it has been fully implemented in the SURFEX platform, specifically as a snowpack scheme within the land surface model ISBA (Noilhan and Planton, 1989;Noilhan and Mahfouf, 1996). This implementation aims particularly at overcoming the limitations mentioned before. The Crocus snowpack scheme is now fully coupled to the ISBA land surface model, allowing straightforward thermodynamic coupling of the snowpack scheme to the soil component of the land surface model. The snowpack scheme Crocus benefits also from coupling routines to several global or regional atmospheric models (GCM: ARPEGE; mesoscale: MESO-NH; mesoscale operational NWP: AROME) as well as facilitated handling of driving data when offline simulations are carried out, including distributed simulations over complex topography. Finally, the implementation of snowpack schemes of varying complexity (e.g. D95, ES and Crocus) within the same land surface model fosters exchanges between model developers and leads to improved capabilities of all models when shared subroutines are improved, thereby minimizing duplication of research work and coding, the latter being prone to errors.
Because several (largely unpublished) evolutions of the scientific content of Crocus have been carried out since its original publications (Brun et al., , 1992, and because the code structure of Crocus in SURFEX has entirely been revisited, this article describes in detail the physical basis and the parameterizations currently implemented in the snow-pack scheme Crocus. It is anticipated that the snowpack scheme Crocus as described here will supersede and replace all previous versions of Crocus developed so far. Our paper is organized as follows. Section 2 gives an overview of the physical processes and the variables included in Crocus. The detailed architecture of Crocus and the physical parameterizations used in the snow scheme are presented in Sect. 3. Section 4 provides technical aspects regarding the format of model inputs and outputs. Validation at a point scale and distributed applications are finally described in Sect. 5.

Physical processes and snow layering
Crocus is a one-dimensional multilayer physical snow scheme. It simulates the evolution of the snow cover as a function of energy and mass-transfer between the snowpack and the atmosphere (radiative balance, turbulent heat and moisture fluxes, ...), and the snowpack and the ground below (ground heat flux). Figure 1 gives an overview of the main physical processes accounted for in Crocus.
The snowpack is vertically discretized on a one dimensional finite-element grid. By convention, the snow layers are described starting from the top of the snowpack to the bottom; the layer number 1 thus corresponds to the surface snow layer (Fig. 1). The vertical discretization is governed by a set of rules, which are designed to develop a realistic dynamic of snowpack layering. These rules are described in Sect. 3.2.
Crocus handles the snowpack stratified parallel to the local slope. The slope angle, referred to as in what follows, has an impact on the compaction rate, since only the component of the weight perpendicular to the snow layering needs to be taken into account. The slope angle also influences the energy and mass fluxes at the snowpack boundaries. As a convention, only vertical incoming and outgoing fluxes are provided to and from the model; the correction of these terms according to the local slope is carried out within SURFEX. Similarly, variables such as total snow depth, total snow water equivalent, and the corresponding variables for each layer are output by the model in terms of their vertical component, i.e. projected vertically.

State variables
In the snowpack scheme Crocus, each snow layer is described by its thickness, D, heat content, H (or enthalpy), density, ρ, and age, A (Fig. 1). Additional variables are used to describe the evolution of snow grains using metamorphism laws (Brun et al., 1992). These variables are dendricity, d, sphericity, s, and grain size, g s . Dendricity characterizes freshly fallen snow and varies from 0 to 1; it roughly represents the remaining initial geometry of snow crystals in the layer, and generally decreases over time in a given layer. Sphericity varies between 0 and 1 and describes the ratio of rounded versus angular shapes. Both variables can be deduced from 2-D image analysis (Lesaffre et al., 1998;Bartlett et al., 2008). An additional historical variable (h) indicates whether there once was liquid water or faceted crystals in the layer. The variables d, s, g s and h are termed the grain variables, and are used to diagnose the snow type (Brun et al., 1992) (Fig. 1). The heat content, H , is used to diagnose the temperature, T , of the snow layer and its liquid water content, W liq (Boone and Etchevers, 2001). Appendix A contains a summary of the variables and units used by the model. The equations governing the evolution of each variable are detailed in the following subsections.

Driving variables
Be it run in coupled or offline mode, the snowpack scheme Crocus within SURFEX needs the following input to run: (i) air temperature, specific humidity and wind speed at a known height above ground; (ii) incoming radiation: direct and diffuse short-wave and long-wave; (iii) precipitation rate, split between rain and snow; (iv) atmospheric pressure. The input for Crocus may be derived directly from local observations, atmospheric models or reanalyses. Section 5 describes several applications of Crocus using different kinds of atmospheric forcing.

Architecture of the snowpack scheme
We only detail here the functioning of the Crocus scheme within SURFEX. Details about SURFEX are provided in Le Moigne et al. (2009). The snowpack scheme Crocus is implemented in SURFEX based on the architecture of the ES snowpack scheme (Boone and Etchevers, 2001). This allows to share common coupling routines between the two schemes. The two main differences between Crocus and ES pertain to the treatment of the vertical grid and the explicit description of snow metamorphism. Other differences regard the parameterizations of physical laws, but the overall structure of the code is similar, as well as the numerical methods used to solve the snow surface/atmosphere exchanges and the set of equations describing the vertical profile of the physical properties of snow. Figure 2 shows an overview of the different calculations performed in the code. Details concerning each process considered by the snowpack scheme Crocus are given in the following subsections, along with the name of the subroutine in charge of the calculations. The routines are described in order of appearance in the code, which corresponds to the chronological order of the computations. Routines which are entirely similar to the ES scheme (Boone and Etchevers, 2001) are not described in detail.

Snowfall
New snowfall is handled by the subroutine SNOWCROFALL UPGRID. When snow is falling, fresh snow layers are added to the snowpack. The model accounts for the impact of near surface meteorological conditions on the properties of falling snow. The density of freshly fallen snow is expressed as a function of wind speed, U , and air temperature, T a , as Geosci. Model Dev., 5, 773-791, 2012 www.geosci-model-dev.net/5/773/2012/ Table 1. Empirical laws for dry snow metamorphism. G is the vertical temperature gradient (|δT /δz|), T the temperature (K) and t is time expressed in days. f , g, h and are empirical functions to predict depth-hoar growth-rate from Marbouty (1980) and are described in Appendix B. where T fus is the temperature of the melting point for water, a ρ = 109 kg m −3 , b ρ = 6 kg m −3 K −1 and c ρ = 26 kg m −7/2 s −1/2 . The minimum snow density is 50 kg m −3 . This density value is then used to convert precipitation amount into snowfall thickness. Variations of density with air temperature and wind speed is plotted in Fig. 3. Parameters in Eq. (1) originate from a study carried out by Pahaut (1976) at Col de Porte (1325 m altitude, French Alps).

Non-dendritic snow
Under strong wind conditions, snowflakes break upon collision between each other and with the snow surface (Sato Table 2. Empirical laws for wet snow metamorphism. θ is the mass liquid water content and t is time expressed in days. v refers to the equivalent volume of snow grain and v ′ 0 and v ′ 1 are empirical constants taken from Brun (1989). Note that θ can be computed from the prognostic variables of the snowpack scheme Crocus as , 2008) so that their properties differ from purely fresh snow (characterized by d = 1 and s = 0.5). Dendricity tends to decrease while sphericity increases. To account for this grain evolution, Guyomarc'h and Merindol (1998) introduced a parameterization which provides dendricity and sphericity of falling snow grains as a function of wind speed, U , (in m s −1 ): Figure 3 presents the dendricity and sphericity of freshly fallen snow as a function of wind speed.
The temperature, hence the heat content of freshly fallen snow, corresponds to snow surface temperature. If no snow is already present on the ground, fallen snow is assigned the minimum value between the ground surface temperature and the temperature of the melting point for water.

Evolution of the vertical discretization of the finite-element grid
The dynamical evolution of the number and thicknesses of the numerical snow layers is a key and original feature of Crocus, which aims at simulating the vertical layering of natural snowpacks in the best possible way (Brun et al., 1992). This feature has been ported into the SUR-FEX implementation of the snowpack scheme Crocus, and is handled by the subroutines SNOWCROFALL UPGRID and SNOWCROGRIDFRESH. The maximum number of numerical layers is an important user-defined set-up option. A minimum number of 3 layers (N min ) is required for solving the heat conduction through the snowpack but there is no limitation on the maximum number. As the maximum number of layers increases, the snowpack stratigraphy can be simulated in more detail. According to the research or operational objectives, the user has to find the appropriate balance between the realism and the computational cost of the simulation. An important point to mention is that the snowpack scheme dynamically manages a different vertical grid mesh, in terms of the number and the thickness of snow layers, for each grid point when it is run in parallel mode for a spatially distributed simulation; this is a common case for snow/atmosphere coupled simulations or for distributed stand-alone simulations. The adjustment of the snowpack layering is achieved with a set of rules within the routine SNOWCROFALL UPGRID. The procedure is activated at the beginning of each time step according to the following sequence: -for snowfall over a bare soil, the snowpack is built up from identical layers, in terms of thickness and state variables. Their number, N, depends on the amount of fresh snow, D new , and on the maximum number of layers, N max : where ⌊.⌋ designates the floor operator.
-for snowfall over an existing snowpack, it is first attempted to incorporate the freshly fallen snow into the existing top layer, provided its grain characteristics are similar and its thickness is smaller than a fixed limit. The similarity between two adjacent layers is determined from the value of the sum of their differences in terms of d, s and g s , each weighted with an appropriate coefficient. If the merging is not possible, a new numerical layer is added to the preexisting layers. If the number of layers then reaches its maximum, a search is carried out to identify two adjacent layers to be merged. This is done by minimizing a criterion balancing the similarity between their respective grain characteristics and their thicknesses; -for no snowfall, a check is carried out to see whether it is convenient to merge too thin snow layers or to split those which are thick. This is achieved by comparing the present thickness profile to an idealized profile, which acts as an attractor for the vertical grid. This idealized thickness profile depends on the current snow depth and on the user-defined maximal number of layers. Figure 4 shows two examples of such an idealized profile. Merging two layers is only possible for those which are similar enough in terms of grain characteristics. Grid resizing affects only one layer per time step, with a priority given to the surface and bottom layers, in order to accurately solve the energy exchanges at the surface and at the snow/soil interface; -for most time steps, no grid resizing is carried out, except that the thickness of each layer decreases according to its compaction rate.
The routine SNOWCROGRIDFRESH ensures the consistency of the physical prognostic variables in case of grid resizing. A projection is achieved from the former vertical grid to the new one. Mass and heat content are conserved. When a new numerical snow layer is built from several former layers, its grain characteristics are calculated in order to conserve the averaged weighted optical grain size of the former layers. This ensures a strong consistency in the evolution of surface albedo, even when frequent grid resizing occurs at the surface in case of frequent snowfalls or surface melting events. Note that the computation of the optical grain size from the snow grain characteristics is detailed in Sect. 3.6.

Snow metamorphism
Snow metamorphism is implemented in a phenomenological way in the snowpack scheme Crocus through a set of quantitative laws describing the evolution rate of the type and size of the snow grains in each layer (Brun et al., 1992). This is carried out within the subroutine SNOWCROMETAMO. A distinction is made between dendritic and non-dendritic snow. Snow initially falls in the dendritic state with dendricity, d, and sphericity, s, given by Eqs. (2) and (3) and remains dendritic until d reaches 0. Snow then reaches the state of rounded crystals, faceted crystals or belongs to an intermediate state. It is then characterized by its sphericity (s), ranging from 0 to 1, and a grain size, g s , ranging from 0.3 to 0.4 mm. Such snow is defined as non-dendritic. The metamorphism laws that govern the evolution of snow grain are given in Tables 1 and 2, respectively, for dry and wet metamorphism. They are similar to the laws initially described by Brun et al. (1992) and are mostly based on empirical fits to experimental data.
Metamorphism laws are used to account for the effect of snow grain type on several parameterizations used to simulate physical process within the snowpack, such as albedo (see Sect. 3.6)

Compaction
The snow layers settle upon the combined effect of snow metamorphism and the weight of the upper layers. The handling of snow compaction is carried out in the subroutine SNOWCROCOMPACTN. The settling is expressed as where D is the layer thickness, σ the vertical stress (computed as the weight of the overlying layers), dt the model time step and η the snow viscosity. The vertical stress from the weight of the overlying layers is computed as follows, for each layer i: where is the local slope, and g is the terrestrial gravitational constant (9.80665 m s −2 ). Note that the vertical stress applied to the uppermost snow layer is equal to half of its own weight. η is described as a function of snow density, temperature, liquid water content, and grain type and is given as follows: where η 0 = 7.62237 10 6 kg s −1 m −1 , a η = 0.1 K −1 , b η = 0.023 m 3 kg −1 and c η = 250 kg m −3 . f 1 and f 2 are two correction factors that adjust the snow viscosity based on snow microstructure properties. They account, respectively, for the decrease of viscosity in presence of liquid water and the increase of viscosity with angular grains: where W liq is the snow layer water content (kg m −2 ), D the snow layer thickness and ρ w the liquid water density, and where g 1 = 0.4 mm, g 2 = 0.2 mm and g 3 = 0.1 mm. Applying Eq. (9) leads to a reduction of the compaction rate in a depthhoar layer.

Impact of wind drift
The compaction and the metamorphism of the surface layers during wind drift events are taken into account in a simplified way, as described in Brun et al. (1997). These calculations are performed within the subroutine SNOWCRODRIFT.
Positive values of S I indicate that snowdrifting occurs while S I = 0 gives the value of the threshold wind speed for snow transport. During a drift event, blown snow particles in saltation break upon collision with the snow surface (Clifton et al., 2006). This results in packing and fragmentation of where τ is empirically set to 48 h. The pseudo-depth in the snow pack, z i (in m, positive downwards), takes into account previous hardening of snow layers j situated above the current layer i: z i = j (D j ×(3.25−S Ij )). Therefore, through the constant Γ drift , compaction and fragmentation rates in a snow layer depend on the grain driftability and are propagated to the layers below with an exponential decay until it reaches a non-transportable layer (S I ≤0). Compaction and fragmentation rates are detailed in Table 3. Brun et al. (1997) introduced this parameterization to simulate a realistic evolution of polar snow density. This turned out to be necessary in polar environments to reproduce correctly the snow thermal conductivity and, therefore, the snow temperature profile (Fig. 3 of Brun et al., 1997). In alpine environments, this parameterization is needed to capture satisfactorily the occurrence of blowing snow events and mass fluxes during those events  As an option and in case of snowdrifting, Crocus computes the associated rate of sublimation following the parameterization developed by Gordon et al. (2006). Under this option, Crocus subtracts the corresponding mass from the snowpack surface. Note that, in stand alone mode, Crocus does not handle explicitly wind-induced snow redistribution since grid points are treated independently from each other. Work is currently in progress to develop the coupling between Crocus and the meso-scale atmospheric model Meso-NH (Lafore et al., 1998) to simulate blowing snow events in alpine terrain.

Snow albedo and transmission of solar radiation
Within the subroutine SNOWCRORAD, the snowpack scheme Crocus handles solar radiation in three separate spectral bands ([0.3-0.8], [0.8-1.5] and [1.5-2.8] µm). First of all, the albedo is computed in each band, as a function of the snow properties in the top 3 cm of the snowpack. In the UV and visible range ([0.3-0.8] µm), snow albedo depends mostly on the amount of light absorbing impurities, but also on its microstructure (Warren, 1982). The latter is represented by the optical diameter of snow, d opt , which corresponds to the diameter of a collection of mono-dispersed ice spheres possessing the same hemispherical albedo as the corresponding semi-infinite snow layer. The impact of the deposition of light absorbing impurities is parameterized from the age of snow. In the near-infrared bands, the spectral albedo depends only on the optical diameter of snow. The optical diameter, d opt , of snow is empirically derived from d, s and g s , based on experimental work by Sergent et al. (unpublished): Once the spectral albedo is calculated, in every spectral band the incoming radiation is depleted by its value, and the remaining part penetrates into the snowpack and is gradually absorbed assuming an exponential decay of radiation with increasing snow depth. The solar flux, Q s , at a depth z below the snow surface is expressed as follows: where R sk represents the incoming solar radiation, α k the albedo and β k the absorption coefficient in the spectral band k. In the current version, the incoming shortwave radiation R s is split into three bands using empirical coefficients (0.71, 0.21 and 0.08, respectively, for band [0.3-0.8], [0.8-1.5] and [1.5-2.8] µm). Future developments will allow to include forcing from an atmospheric model where incoming shortwave radiation is partitioned into several bands. Shortwave radiation excess for thin snow cover (transmitted through the snow) is added to the snow/ground heat flux. The albedo and the absorption coefficient for each spectral band are given in Table 4.

Surface fluxes and surface energy balance
The routine SNOWCROEBUD calculates the aerodynamic resistance and the turbulent exchange coefficients between the snow surface and the atmosphere following the same approach as Boone and Etchevers (2001). Those coefficients are then used by SNOWCROFLUX to compute surface fluxes. The latent heat flux includes contributions from evaporation of liquid water in the surface layer and sublimation. It is written as where L f and L v denote the latent heat of fusion and vaporisation, respectively, q a is atmospheric specific humidity Geosci. Model Dev., 5, 773-791, 2012 www.geosci-model-dev.net/5/773/2012/ Table 4. Evolution of snow albedo and absorption coefficient for three spectral bands based on theoretical studies of Warren (1982). A is snow-surface age expressed in days and d opt (m) the optical grain diameter given by Eq. (13). The term P /P CDP represents the decreasing effect of ageing on the albedo with elevation (P : mean pressure and P CDP : 870 hPa). (kg kg −1 ), q sat (T ) is the saturation specific humidity above a flat ice surface at the temperature T and T s is snow surface temperature. χ denotes the ratio between the solid and liquid phases of the turbulent mass exchanges between the snow surface and the atmosphere. It is evaluated in SNOWCROEBUD, according to the following rule: the absence of liquid water in the surface layer at the beginning of the time step imposes only solid exchanges (hoar deposition or sublimation); the presence of liquid water imposes liquid condensation or evaporation; in the case where the computed evaporation leads to the complete removal of the available liquid water, the ratio between the solid and liquid phases is adjusted in order to extract the remaining mass from the ice. The sensible heat flux is where C p is the specific heat of air and s and a are Exner functions for the surface and the atmosphere, respectively. The formulation of the turbulent exchange coefficient C H follows Noilhan and Mahfouf (1996) and is based on Louis (1979): where z u and z a are the heights of the wind and air temperature measurements and κ is the von Karman constant. The effective roughness z 0 takes into account the effects of both snow and vegetation. f (R i ) represents the dependence of the transfer coefficient on the atmospheric stability (function of the Richardson number, R i ). In contrast to the first version of Crocus (Brun et al., 1992;Martin and Lejeune, 1998), C H is not treated as a site-specific calibration parameter. However, as Martin and Lejeune (1998) suggest, C H values can, under certain conditions, still become quite low, thereby effectively decoupling (too much) the surface from the atmosphere. A model option exists which consists of the use of a maximum Richardson number (R i max ) for very stable conditions. The incorporation of an effective roughness z 0 is especially important for local studies near or within forest or in a spatially distributed simulation with vegetated areas within the computational cells. ISBA partitions the grid cell between vegetation and bare ground. Both of them may be covered by snow with expressions of fractional snow covered area (FSCA) calculated from SWE and vegetation roughness (Douville et al., 1995). FSCA is then used to compute the effective roughness and to partition the flux of heat, momentum and mass between the snow and non-snow covered fractions of the grid cell. Distributed applications of the model require such an approach in order to represent snow cover heterogeneity within a grid cell. However, for point scale applications focusing on snow physics, an option in SURFEX forces FSCA to 1 as soon as the snowpack reaches a relatively low user-defined SWE threshold. This option is recommended for local scale applications with an emphasis on studying snow physics such as the simulations carried out at Col de Porte (see Sect. 5.1).

Resolution of snow temperature profile
The heat diffusion within the snow cover is computed by SNOWCROSOLVT using the implicit backward-difference integration scheme of ISBA-ES (Boone and Etchevers, 2001). The snow effective thermal conductivity, k, expressed in W m −1 K −1 follows the expression of Yen (1981): The net heat flux at the snow-atmosphere interface combines the turbulent fluxes (described in the previous section) with the net radiative components (short-and longwave). It also includes a precipitation heat advection term when it is raining for offline local studies. In terms of longwave radiation, the snow emissivity is assumed to be 1. that is generally between one to several centimeters thick, depending on the local soil characteristics and on the soil scheme options. We recommend to use the version of ISBA based on a multi-layer diffusive approach [ISBA-DF] (Boone et al., 2000;Decharme et al., 2011) to simulate the evolution of the soil temperature and water content (both liquid and ice). ISBA requires the knowledge of the soil texture (fractions of root, clay, sand and silt). They can be provided by the user for point specific simulations or taken from global database available at 1-km resolution for distributed simulations (ECOCLIMAP, Masson et al., 2003). The flux calculation differs from the first version of Crocus (Brun et al., , 1992 where the ground heat flux was imposed depending on the geographic location, the elevation, and the season. When coupled to an atmospheric model, SURFEX has a user option to couple the snow and the atmosphere using a fully implicit numerical coupling (Polcher et al., 1998;Best et al., 2004). It is especially adapted for relatively large time steps, such as those used for long range NWP or GCM experiments. The model can be run using time steps up to 1 hour in offline mode and 30 minutes when coupled to a GCM.

Snow melt
Total or partial melting of the snowpack is handled by three subroutines: SNOWCROGONE, SNOWCROLAYERGONE and SNOWCROMELT.
SNOWCROGONE inherits ISBA-ES features. It calculates the new heat content of the snowpack from the new temperature and density profile. It compares this energy to the amount of energy which is necessary for the complete melting of the snowpack ice mass, from which possible sublimation has been subtracted. If the available energy exceeds this energy, the snowpack completely melts and the routine computes the corresponding impact on the ground heat and water fluxes, in order to ensure the conservation of energy and mass, while taking into account the vapor exchanges between the vanishing snowpack and the atmosphere.
SNOWCROLAYERGONE accounts for the case when one or several snow layers completely melt during a time step, before the computation of the partial melting/refreezing inside each snow layer. First, the routine compares the new heat content of each snow layer to the amount of energy which is necessary for the complete melting of its ice mass. Then, if the available energy exceeds it, the snow layer is merged with the underlying layer, except for the bottom layer which is merged to the overlying layer. Each new merged layer conserves the energy and mass of the two layers it is made from. It inherits the grain size, shape, history and age from the layer with which the melted layer has been merged.
SNOWCROMELT is run after SNOWCROGONE and SNOWCROLAYERGONE, which means that the available energy from the new temperature of any snow layer is not large enough to melt it completely. Then, when the new temperature of a layer exceeds the melting point, the tem-perature is turned to the melting point and the corresponding energy is consumed for ice melting. The corresponding melt water is added to the liquid water content of the layer. The dry density of melting layers is conserved at this stage and their thickness decreases accordingly.

Water flow and refreezing
The routine SNOWCROREFRZ handles the refreezing of liquid water and its flow through the snow pack. It first updates the liquid water content of the surface snow layer by including contributions from rainfall and liquid condensation or evaporation at the surface (calculated by SNOWCROEBUD). Then, it calculates the amount of energy available for liquid water freezing from the new temperature of each snow layer. If freezing occurs in a given layer, its liquid water content is decreased and its temperature is modified accordingly. The water flow through the snow layers is then simulated. The liquid water content of the snowpack is modeled as a series of reservoirs (one for each layer). Water flow occurs when the liquid water content exceeds the maximum liquid water holding capacity (W liq max in kg m −2 ). It is expressed as 5 % of the total pore volume (Pahaut, 1976): where ρ w and ρ i are the water and ice density, respectively. The model considers only gravitational flow and neglects the formation of capillary barriers (Jordan, 1995). Water leaving the bottom of the snowpack is available to the soil for infiltration or surface runoff. The water flow solution procedure starts from the uppermost layer and proceeds downward. Water entering a layer refreezes if thermodynamics allows it. Once a layer can no longer freeze liquid water present in the layer (i.e. T = T fus ), then the unfrozen water is retained up to the maximum holding capacity. The refreezing and water retention processes increase the layer-average density and mass. Water flow processes do not impact the layer thicknesses.

Snow sublimation and hoar deposition
The routine SNOWCROEVAPN adds or substracts to the snow surface layer the ice amount corresponding to the turbulent vapor fluxes, according to the ratio between the solid and liquid phases which have been determined in SNOWCROEBUD. The surface snow layer thickness is adjusted accordingly while the density is assumed to stay unchanged. This implies that at this stage of development, the snowpack scheme Crocus does not represent the specific properties of surface hoar.

Final updates: surface albedo, heat content
The final updates ensure the coherence between the final snowpack properties and the variables stored at each time step: Geosci. Model Dev., 5, 773-791, 2012 www.geosci-model-dev.net/5/773/2012/  -the new snow albedo depends on the final snow grain type close to the surface following Table 4, -the heat content H for each layer is computed using the current snow temperature and liquid water content.

Format of model input/output
This section deals with model input/output in the context of forced simulations (offline). Indeed, coupled simulations are driven by the atmospheric model which generally also handles the model output.
Except formats specific to atmospheric models, the main formats for driving data are ASCII, binary and NetCDF (Rew and Davis, 1990). The latter is preferred for distributed simulations over many points, as its data structure is dedicated to handling multi-dimensional datasets easily (Zender, 2008).
Model output can be provided in various formats, but we only describe here the (recommended) use of the NetCDF output. Model output settings are a general feature of SUR-FEX, thus there is no dedicated model output in the case of snowpack simulations. Data relevant to the snowpack state are provided in two output files at the level of the ISBA land surface scheme within SURFEX. The first one, termed ISBA PROGNOSTIC.nc contains the values of the state variables of the snow and ground layers. The second one, termed ISBA DIAGNOSTIC.nc, contains diagnosed quantities such as surface fluxes, albedo, surface temperature, melt water runoff etc. The main dimension of both files is the time and the location. The latter can either be two-dimensional (rectangular regular grid, lat/lon or x/y) or one-dimensional. Specific routines are used after a model run using the snowpack scheme Crocus, adding the dimension snow layer to the prognostic output file, i.e. for each time and location, each snow variable is then represented as a single data vector. The resulting data file in NetCDF follows an ad-hoc, hitherto internal format termed the "NetCDF Snowpack Profile Format". This data format aims at complying with the NetCDF Climate and Forecast (CF) Metadata Convention (Gregory, 2003). However, because the thickness and the number of snow layers vary in time, there is no fixed vertical grid for storing the vertical profiles. Instead, the thickness of each layer is provided as a fully fledged output variable: the data SNOWDZ (snow layer thickness, in m) in such a file has dimension (time, snow layer, location), where the dimension snow layer starts from the uppermost snow layer downwards and contains the maximum number of snow layers considered in the simulation (in general, 20 or 50). In the case where the maximum number of snow layers is not reached, "empty" layers are treated as missing data using the NetCDF standard practice. Other variables (snow temperature, liquid water content, etc.) are stored accordingly. Data in this file can also be vertically integrated for variables, such as snow depth, SWE, or uppermost soil layer temperature or liquid water content. The use of this data format greatly facilitates data storage, handling, post-processing (including plotting) and further computations from the model output, such as mechanical stability evaluation using, e.g. the MEPRA algorithm (Durand et al., 1999) or coupling to microwave emission models (Brucker et al., 2011). Dedicated tools for the plotting of individual snowpack profiles or temporal overviews of the time evolution of the physical properties of snow are being developed from this data format. An example of the time evolution of the internal physical properties of snow is provided in Fig. 5, based on one year of model output from the model run at Col de Porte, France, described in detail in Sect. 5.1.

Model evaluation and examples of use
The following sections present model runs used to evaluate Crocus within SURFEX, as well as providing illustrations of the versatile use of this new implementation of this snowpack scheme. Note, however, that the development of the snowpack scheme Crocus within SURFEX benefited from earlier experience with both the ES snowpack scheme (Boone and Etchevers, 2001) and the Crocus snowpack model (Brun et al., , 1992 architecture, so that no large difference in model behavior was anticipated. Nevertheless, the examples shown below demonstrate that the snowpack scheme Crocus within SUR-FEX behaves quite similarly, and generally better than the original Crocus snowpack model.

Offline simulation and detailed evaluation of 18 snow seasons at the Col de Porte site (1993-2011)
The meteorological station at Col de Porte (1325 m altitude, 45 • 17 ′ N, 05 • 45 ′ E) in the Chartreuse mountain range near Grenoble, France, has been used for over 50 years as an experimental field site devoted to the study of snow in mountains. Data for driving and evaluating snowpack models have been collected at the appropriate time scales for several decades. Data from the Col de Porte (CDP) have thus been widely used in the past for model development and evaluation, such as the original Crocus snowpack model (Brun et al., , 1992 and the international Snow Model Intercomparaison Project (SnowMIP) initiative (Etchevers et al., 2004). We here present a single model run carried out using the snowpack scheme Crocus within SURFEX, using driving data from CDP. Much of the focus of studies carried out at CDP is on the snow season, thus meteorological data are quality-controlled for the periods of time when snowfall happens, i.e. from 20 September to 10 June of the following year. To perform a continuous run without data-gap during the summer, we use the output of the SAFRAN downscaling tool to provide meteorological driving data to the land surface model from 10 June to 20 September of each year (Durand et al., 1999). Using quality-controlled data from the CDP in-situ meteorological data for the snow season, a single forcing data file was built, covering the period between 1 August 1993 to 31 July 2011. It consists of hourly records of the driving data for the land surface model ISBA within SURFEX. Full details regarding the dataset and its availability are given in Morin et al. (2012a). The model run was initialized with no snow on the ground on 1 August 1993, and a single run was performed until 31 July 2011. The soil configuration corresponds to the multilayer diffusion scheme (ISBA-DF) (Decharme et al., 2011), where 20 soil layers were considered down to a depth of 10 m below the surface. The run using the Crocus snowpack scheme was carried out allowing up to 50 snow layers. A similar model run was carried out using the intermediate complexity snow scheme ES (Boone and Etchevers, 2001) instead of Crocus. In addition, model runs were performed using the same driving data and the original snowpack model Crocus as described in Brun et al. (1992). In the case of the snowpack schemes coupled to ISBA (ISBA-Crocus and ISBA-ES), the model runs were carried out by setting a snow fraction of 1. (see Sect. 3.7 for details), and an effective roughness length z 0 = 5 mm. This value corresponds to a near-optimum for both models, which can be viewed as a consequence of the fact that they share a similar surface energy budget, although the physics within the snowpack are different (more detailed in ISBA-Crocus). Both model runs were evaluated against daily observations of snow depth and  (Kodama et al., 1979;Paquet and Laval, 2006), with an uncertainty on the order of 10 %. For consistency reasons between the two records of evaluation data, the simulation was evaluated concurrently for snow depth and SWE for the ten winter seasons between 2001 and 2011. Figure 6 shows an overview of simulations in terms of snow depth and SWE and the corresponding statistics in terms of bias and root-mean-square-deviation (rmsd). Table 5 provide statistics computed using the full 10 yr period of simulation. A previous intercomparison between Crocus and ISBA-ES at CDP indicated that SWE and snow depth were significantly better simulated by ISBA-ES than Crocus at this site for one season (see e.g. Table 2 of Boone and Etchevers, 2001). The same statistics computed for 10 snow seasons at CDP indicate large year-to-year variations in model performance, preventing a fully informative comparison between model skills based on a limited amount of model years. Our interpretation of the statistics computed is that all three models (Crocus, ISBA-Crocus and ISBA-ES) perform satisfactorily at CDP in terms of bulk snowpack properties. A more detailed investigation of ISBA-Crocus performance in terms of physical properties of snow (density and microstructure) was recently carried out by Morin et al. (2012b).
The simulation also provides the information that the computational cost of running the snowpack scheme Crocus with a maximum of 50 numerical snow layers is only 2.3 times larger than for running the snowpack scheme ES, which remains on the same order of magnitude as for previous such comparisons (Boone and Etchevers, 2001, factor 2.5).

Distributed offline simulation of the snowpack at the spatial scale of a mountain range in the Alps
The general framework of SURFEX permits spatially distributed simulations over a given domain. Here we present the example of the evolution of the snowpack over the Grandes Rousses moutain range in the French Alps during . The distributed simulation is based on a digital elevation model with an horizontal resolution of 150 m, which allows a fine representation of the differences in terms of radiative budget between the simulation points.
The meteorological forcing was based on the output from SAFRAN (Durand et al., 1999) over the Grandes Rousses range, i.e. hourly meteorological driving data (Sect. 2.3) for six different aspects (N, E, SE, S, SW, W) at 300 m elevation intervals. This information was interpolated to each grid point as a function of its elevation, local slope and aspect (Fig. 7). Incoming shortwave radiation was corrected to account for effects of slope aspect and terrain shading.
The simulation started from 1 August 2010 over a snowfree domain, and lasted until 1 May 2011. Wind-induced snow transport was not explicitly included. Figure 7 (3) shows a map of snow depth over the simulation domain. Strong contrasts are observed in terms of snow depth between the north-facing and south-facing slopes due to topographic effects on the surface energy balance.

Atmosphere/snow coupled simulation of the energy balance of the snowpack in Antarctica
One of the first applications of the implementation of Crocus into SURFEX has been the set-up and the evaluation of a 11-day detailed 3-D coupled snow/atmosphere simulation around Dome C (Brun et al., 2011). From a technical point of view, the set-up of such a configuration has been extremely simplified by the general SURFEX framework, which includes the algorithms and interfaces allowing a full-coupling between the different surface schemes and the atmosphere. It was based on a configuration of the AROME regional meteorological model , over a 625 × 625 km 2 domain centered around Dome C, Antarctica. The horizontal resolution was 2.5 km and 60 vertical levels were used, allowing a very detailed vertical resolution in the lower layers of the atmosphere. The snow model included 20 snow layers, representing the top 10-m of the firn and snowpack, initialized from local observations. The evaluation was based on a comparison between the observed and simulated snow temperature profiles, and temperature and wind profiles in the atmospheric boundary layer. In spite of a poor simulation at times of clouds, the surface and near-surface snow temperatures were correctly simulated (Figs. 8 and 9), showing neither significant bias nor drifts during the simulation period. This study proved to be very encouraging for improving the detailed representation of the physical processes at the snow/atmosphere interface, either in climate models or in NWP systems.

Conclusions
This paper describes the new version of the snowpack scheme Crocus. It includes the main features of the previous versions of Crocus in terms of dynamical layering of the snowpack and explicit representation of snow metamorphism (Brun et al., , 1992. The surface energy balance and heat redistribution within the snowpack are now solved following the ISBA-Explicit Snow (ES) snowpack scheme (Boone and Etchevers, 2001). New parameterizations such as the impact of wind-drift allow Crocus to be run in different environments from polar regions to alpine terrain. This version of Crocus has then been implemented within the surface module SURFEX to better represent the interactions between the snowpack and its environment. Crocus is indeed fully coupled to the ISBA land surface model and its soil component allowing for an accounting of snow-vegetation interactions in a simplistic manner and realistic soil heat flux below the snow cover. As a general platform used by Obs -43cm Sim -43cm Obs -63cm Sim -63cm Obs -83cm Sim -83cm Météo France NWP and climate models, SURFEX can be coupled to several atmospheric models. Therefore, the snowpack scheme Crocus can be run either in stand alone mode, using a time series of meteorological forcing (single point or distributed), or in a fully-coupled mode (explicit or fully implicit) with an interactive simulation of the atmosphere. This enables Crocus to be used for many applications including avalanche forecasting, hydrological or climate studies. A 10-yr evaluation (2001-2011) of the new snow scheme has been carried out at the Col de Porte experimental site (French Alps). Results show that ES and Crocus perform well and with comparable levels of performance, in terms of snow depth, SWE and numerical costs. When coupled to the atmospheric model AROME over Dome C (Antarctica), Crocus was able to reproduce reasonably well the evolution of the snow surface temperature over an 11-day period (Brun et al., 2011). The coupling of the atmospheric model with Crocus/SURFEX also proved to be able to simulate a consistent evolution of the atmospheric boundary layer. In alpine terrain, model applications include the simulation of the sea-sonal evolution of the snowpack over a whole mountain range using distributed meteorological forcing.
Further developments of the snowpack scheme Crocus within SURFEX are planned. In terms of the snowpack scheme itself, the two major planned developments are the comprehensive revisit of the solar radiation transfer scheme, and the reformulation of the snow metamorphism laws. Beyond the scope of the snowpack scheme Crocus itself, an explicit representation of snow/canopy interactions is currently being developed within ISBA. This will permit an explicit representation of turbulent and radiative transfer within and below the canopy, and certain processes critical for modeling snow in a forest, such as unloading. The coupling of Crocus with the atmospheric model Meso-NH is also in progress and will lead to the inception of a modeling platform dedicated to the simulation of the snowpack evolution during snow-drift events.
The snowpack scheme Crocus fully belongs to SURFEX from version 7 on, and is available for research purposes on request to the authors. Information pertaining to the evolution