Simulation of wind-induced snow transport and sublimation in alpine terrain using a fully coupled snowpack/atmosphere model

In alpine regions, wind-induced snow transport strongly influences the spatio-temporal evolution of the snow cover throughout the winter season. To gain understanding on the complex processes that drive the redistribution of snow, a new numerical model is developed. It directly couples the detailed snowpack model Crocus with the atmospheric model Meso-NH. Meso-NH/Crocus simulates snow transport in saltation and in turbulent suspension and includes the sublimation of suspended snow particles. The coupled model is evaluated against data collected around the experimental site of Col du Lac Blanc (2720 ma .s.l., French Alps). First, 1-D simulations show that a detailed representation of the first metres of the atmosphere is required to reproduce strong gradients of blowing snow concentration and compute mass exchange between the snowpack and the atmosphere. Secondly, 3-D simulations of a blowing snow event without concurrent snowfall have been carried out. Results show that the model captures the main structures of atmospheric flow in alpine terrain. However, at 50 m grid spacing, the model reproduces only the patterns of snow erosion and deposition at the ridge scale and misses smaller scale patterns observed by terrestrial laser scanning. When activated, the sublimation of suspended snow particles causes a reduction of deposited snow mass of 5.3 % over the calculation domain. Total sublimation (surface + blowing snow) is three times higher than surface sublimation in a simulation neglecting blowing snow sublimation.


Introduction
Wind-induced snow transport is an important component of the interaction between the cryosphere and the atmosphere. It occurs in regions seasonally or permanently covered by snow. In alpine terrain, snow transport creates inhomogeneous snow depth distribution, strongly influenced by the local topography (e.g. Durand et al., 2005;. Snow is eroded in areas exposed to strong wind (crest for example) and is deposited in areas sheltered from the wind. This has a major influence on the evolution of the avalanche danger (Schweizer et al., 2003). Cornices and wind slabs result indeed from the deposition of small rounded grains during blowing snow events and their formation affects the snowpack stability (Meister, 1989). Windinduced snow transport has also hydrological consequences. Part of the transported snow mass is indeed lost by sublimation (e.g. Groot Zwaaftink et al., 2011, hereafter referred to as GZ11) and differential melt patterns are produced by drifts and scours affecting the springtime runoff patterns (e.g. Winstral et al., 2013).
Wind-induced snow transport occurs when the wind speed exceeds a threshold value that depends on the snow type at the surface (e.g. Guyomarc'h and Mérindol, 1998). Three modes of snow transport are generally identified: reptation, saltation and turbulent suspension. Reptation corresponds to the rolling of particles over the surface of the snowpack and its contribution is negligible compared to the other processes (Kosugi et al., 1992). It is commonly neglected in blowing snow models. In saltation, particles follow ballistic trajectories in a shallow layer close to the ground. When returning to the surface, they may rebound and/or eject new grains. Turbulent suspension occurs in the atmosphere above the saltation layer where snow grains are transported by turbulent eddies without contact with the surface. Distances of transport are limited by the sedimentation and sublimation of snow grains. The latter process modifies the vertical profiles of temperature and humidity in the surface boundary layer (Schmidt, 1982;Déry et al., 1998).
Several models have been developed to simulate windinduced snow redistribution in alpine terrain and resulting snow-depth pattern (e.g. Naaim et al., 1998;Gauer, 1999;Durand et al., 2005;Liston et al., 2007;Lehning et al., 2008;Schneiderbauer and Prokop, 2011). These models are generally made of two components: (i) a snowpack component to estimate the threshold wind speed for snow transport and the erodible snow mass and (ii) an atmospheric component to simulate at high resolution the spatial and temporal evolution of the wind field and the resulting snow transport. They cover a wide range of complexity depending if they are dedicated to simulate single blowing snow event or the entire snow season (Liston et al., 2007).
Simulating wind-induced snow transport in alpine terrain requires primarily a good knowledge of the high-resolution wind field over complex topography. Crest speed-up, flow channelling or the formation of recirculation zone are indeed the driving mechanisms behind the inhomogeneous snow distribution resulting from blowing snow events (e.g. . To reproduce such features, the more advanced snow transport models use three-dimensional wind fields. They are computed by models of computational fluid dynamics (Gauer, 1999;Schneiderbauer and Prokop, 2011) or atmospheric models run in the Large Eddy Simulation (LES) mode. Bernhardt et al. (2010) combined for example wind fields from the MM5 atmospheric model (Grell et al., 1995) at an horizontal resolution of 200 m with the snowtransport model of Liston et al. (2007). They used kinematic downscaling to refine MM5 wind fields from 200 to 30 m. The atmospheric model ARPS (Xue et al., 2000) provided wind fields at an horizontal grid spacing down to 5 m  to drive the snow transport module of Alpine3D .
In alpine terrain, the atmospheric models previously mentioned are used to drive snow transport models but do not simulate drifting and blowing snow in a coupled mode. The driving wind fields are extracted from wind field libraries and kept constant for a given time step (1 h for example in the model of Lehning et al., 2008). Until now, the coupled simulation of snow transport has been only implemented in regional climate models applied over large areas at horizontal resolution higher than 5 km (Gallée et al., 2001;Lenaerts et al., 2012). However, previous works have shown that atmospheric models can be run at high resolution in complex terrain to simulate in coupled mode meteorological situations such as wildland fire (Mandel et al., 2011, resolution of 100 m) or scalar dispersion (Michioka and Chow, 2008, resolution of 25 m). These studies were successful at capturing the flow structures in complex terrain. As a consequence, atmospheric models can be applied to the coupled simulation of blowing snow events in alpine terrain.
In this study, we introduce a new coupled model to simulate blowing snow events in alpine terrain and resulting snow redistribution. The atmospheric component is made of the atmospheric model Meso-NH (Lafore et al., 1998) to simulate the evolution of meteorological conditions and resulting snow transport. At the bottom of the atmosphere, the latest version of the detailed snowpack model Crocus  describes snowpack properties. Examples of problems that can be better investigated with fully coupled snowpack/atmosphere simulations of blowing snow events in alpine terrain include (1) the importance of blowing snow sublimation and its feedback on the atmospheric boundary layer (GZ11), (2) the use of grid nesting techniques to provide realistic boundary conditions to local atmospheric simulation (resolution ≈ 50 m) (e.g. Talbot et al., 2012) and (3) the relative contribution of preferential deposition of snowfall and wind-induced snow transport to the spatial variability of the mountainous snowcover .
Our paper describes the new snow transport model Meso-NH/Crocus. It is evaluated against a set of observations collected during two blowing snow events at the experimental site of Col du Lac Blanc (French Alps). Our paper is organised as follows. Section 2 describes the two components of the coupled snowpack/atmosphere model. Then, Sect. 3 presents the new drifting and blowing snow scheme included in the coupled model. Section 4 gives additional information regarding Col du Lac Blanc and the evaluation strategy. Model results in 1-D configurations are shown in Sect. 5. The 3-D simulation of a blowing snow event is then evaluated in Sect. 6. The influence of blowing snow sublimation is also presented in this section. Finally, model results are discussed in Sect. 7.

Atmospheric model
We use the mesoscale, nonhydrostatic atmospheric model Meso-NH. This model has been jointly developed by CNRM (Météo-France) and Laboratoire d'Aérologie (CNRS) (Lafore et al., 1998). Meso-NH can simulate fine scale (LES type) to synoptic scale (horizontal resolution ranging from a few metres to several tens of kilometres) and can run in two-way nested mode. For high-resolution simulations, the scheme of Cuxart et al. (2000) allows the computation of 3-D turbulence. Cloud dynamics and resulting precipitation (rainfall and snowfall) are simulated by the bulk microphysical scheme of Pinty and Jabouille (1998 Meso-NH has been previously applied in alpine terrain during the Mesoscale Alpine Programme (e.g. Stein, 2004). Studies were carried out at sub-kilometre resolutions (0.3-1 km) and focused on Foehn flow (Beffrey et al., 2004) or orographic precipitations (Lascaux et al., 2006). More recently, Amory (2012) showed that Meso-NH is suitable to simulate high-resolution wind field (down to 12 m) over complex terrain.

Snow model
Meso-NH is coupled to the surface platform SURFEX which handles energy and mass fluxes between the atmosphere and the surface (Masson et al., 2013). SURFEX includes in particular the land surface scheme ISBA (Interactions between Soil, Biosphere, and Atmosphere, Noilhan and Planton, 1989) which represents the snowpack through several schemes of increasing complexity. In this study, we use SURFEX/ISBA-Crocus (referred later as Crocus), the latest version of the snowpack scheme Crocus (Brun et al., 1989(Brun et al., , 1992 which has been recently implemented in SUR-FEX . Crocus has been run operationally for avalanche forecasting over the French mountains for 20 yr (Durand et al., 1999).
Crocus is a one-dimensional multi-layer snowpack scheme. It simulates the evolution of the snowpack as a function of energy and mass-transfer between the snowpack and the atmosphere (radiative balance, turbulent heat and moisture transfer, precipitation) and the snowpack and the ground below (ground heat flux). One important feature of the model is its ability to simulate snow metamorphism through a comprehensive set of semi-empirical laws. The type and size of the crystals of each layer of the snowpack are prognostic variables, which control the surface albedo and the compaction rate of the different layers. A complete description of the model can be found in Vionnet et al. (2012).

Treatment of drifting and blowing snow
We have implemented new specific routines within the coupled system Meso-NH/Crocus to handle wind-induced snow transport. As in earlier-developed models (e.g. Gauer, 1999;Durand et al., 2005;Liston et al., 2007;Lehning et al., 2008), snow transport is divided into turbulent suspension and saltation. Figure 1 shows an overview of the new blowing snow scheme which is made of three main components: (i) a specific scheme implemented in Meso-NH simulates snow transport in turbulent suspension, (ii) dedicated routines in SURFEX determine blowing snow occurrence and simulate snow transport in saltation and (iii) mass exchange between the snowpack and the atmosphere are computed through an adapted version of the surface boundary layer (SBL) scheme Canopy (Masson and Seity, 2009). Details concerning each component are given in the following subsections while Ap-pendix A contains a summary of the variables and units used by the model.

Double-moment scheme
Several field experiments have shown that the size distribution of blown snow particles in the atmosphere is satisfyingly represented by a two-parameter gamma distribution (Budd, 1966;Dover, 1993;Nishimura and Nemoto, 2005;Naaim-Bouvet et al., 2011). The snow particle size distribution (PSD) follows: where N s is the total number of particles per unit of mass (kg −1 ), ρ air the air density (kg m −3 ), r the particle radius (m) and Ŵ the Gamma function. α (-) and β (m) denote the shape and scale parameters of the distribution. The average radius, r m , is given by r m = αβ.
In our model we consider blown snow particles as spherical and use a double-moment scheme to represent the temporal and spatial evolution of the two-parameter gamma distribution of snow in suspension (Déry and Yau, 2001). Such approach predicts the evolution of the number concentration per unit of mass, N s , and mixing ratio, q s , of blown snow particles and requires a specified value of α. At Col du Lac Blanc, Naaim-Bouvet et al. (2011) have reported values of α ranging from 3 to 4. In the rest of the paper, we take α = 3 for the use of Meso-NH/Crocus in alpine terrain.
The equations for N s and q s depend on space variables x i and time t according to: where Adv denotes the advection term and Turb the turbulence term. Sedim indicates the sedimentation term and Subl represents the sink term associated to sublimation for N s and q s . u is the 3-D wind vector, V N and V q are the number-and mass-weighted mean particle fall speeds. The advection of N s and q s is handled by specific routines of Meso-NH dedicated to the advection of meteorological variables. It relies on the Piecewise Parabolic Method (Colella and Woodward, 1984;Carpenter et al., 1990). The flux limiter of Skamarock (2006) ensures monotonicity preservation.
The 3-D turbulence scheme of Meso-NH (Cuxart et al., 2000) computes the turbulent diffusion of N s and q s . This scheme solves a prognostic equation for the turbulent kinetic energy, e TKE . Turbulent fluxes are then computed following a 1.5 order closure. In our model, turbulent fluxes of blowing snow variables are proportional to the turbulent fluxes of scalar variables. We define ζ as the ratio between the diffusion coefficient of scalar variables, K Sca , and blowing snow variables, K Snw . Thus, the turbulent flux q ′ s u ′ j is written: where C S and the stability functions j are taken from Redelsperger and Sommeria (1981). L m denotes the mixing length. For 3-D simulations at high horizontal resolution (≤ 200 m), it depends on the mesh size and is limited by the mixing length of Deardorff (1972) in stable cases. Close to the surface, L m follows the formulation of Redelsperger et al. (2001). A sensitivity study showed that the value ζ = 0.25 allows to reproduce vertical profiles of blowing snow fluxes measured at Col du Lac Blanc over a large range of wind speeds as discussed in Vionnet (2012). Sedimentation terms in Eqs.
(2) and (3) require the computation of the number-weighted and mass-weighted mean particle fall speeds. These velocities depend on the PSD following: where v(r) denotes the fall speed for a particle of radius r. v(r) follows a quadratic equation derived by Dover (1993) based on a balance between the gravitational force acting on a spherical particle and the drag force using the drag coefficient proposed by Carrier (1953): and B = 5.516ρ ice 4ρ air g where ν air (m 2 s −1 ) is the air viscosity and g (m s −2 ) the acceleration due to gravity. This expression for v(r) represents the transition from laminar regime where Stokes Law applies (1 ≤ r ≤ 30 µm) to turbulent regime. It is therefore suitable for blowing snow particles since it covers their typical range of radius (20-200 µm, e.g. Budd, 1966). Using Eq. (7) for v(r), no analytical solution can be found for V N and V q so that the integrals in Eqs. (5) and (6) are solved numerically. To save computational time when 3-D simulations are performed, a model option allows to use pre-computed look-up tables of V N and V q function of the average radius r m and the air pressure P air .

Blowing snow sublimation
Sublimation terms appear in Eqs.
(2) and (3). Indeed, when transported, suspended snow particles undergo sublimation if the ambient air is unsaturated with respect to ice. The sublimation of blowing snow acts as a sink of snow mass, a source of water vapour and a sink of sensible heat in the atmosphere. Several blowing snow models calculate sublimation rates of blown snow particles and account for the loss of mass due to sublimation (e.g. Déry et al., 1998;Bintanja, 2000;Liston et al., 2007, GZ11). They all rely on the formulation of Thorpe and Mason (1966) Wever et al. (2009) showed that this formulation can be transferred to an ensemble of snow particles. Our model follows the same formulation. Ignoring the influence of solar radiation, the sublimation rate of an ice sphere of radius r is: where σ = (e − e si )/e si is the water vapour deficit with respect to ice, with e and e si the vapour pressure and its value at saturation over ice (Pa) at air temperature T air (K).
R v denotes the gas constant for water vapour (J kg −1 K −1 ), D v is the molecular diffusivity of water vapour in air (m 2 s −1 ), K air is the molecular thermal conductivity of the air (J m −1 s −1 K −1 ), L s is the latent heat of sublimation (J kg −1 ). Nu and Sh are, respectively, the Nusselt and Sherwood numbers. Expressions for D v and K air are taken from Pruppacher et al. (1998). The Nusselt and Sherwood numbers represent the heat and mass transfer between the ice particle and the atmosphere. They are related to the particle Reynolds number Re (Lee, 1975): where V v is the ventilation speed defined as the relative speed between the snow particle and the air. It is taken equal to the settling velocity of the particle given by Eq. (7) (Schmidt, 1982;Déry et al., 1998). This approach neglects ventilation effects associated with turbulence (Dover, 1993;Bintanja, 2000). The total mass sublimation rate, S q (kg kg −1 s −1 ), in Eq. (3), is then obtained through integration over the particle size spectrum (Déry and Yau, 1999): where Nu(r m ) is computed as in Eq. (9) with r m for the radius and V q for the settling velocity. The total number sublimation rate, S N (kg −1 s −1 ), is given by: S N = S q N s /q s (Déry and Yau, 2001). Sublimation feedbacks on the air are represented through additional terms proportional to S b in the prognostic equations for vapor mixing ratio and potential temperature (Déry and Yau, 2001, GZ11).

Occurrence of snow transport
The occurrence of snow transport is determined by the formulation of Guyomarc'h and Mérindol (1998) which gives the 5 m threshold wind speed, U 5t , as a function of snow grain type at the surface. The presence of a wet layer or a crust layer at the top of the snowpack prevents snow drifting. Vionnet et al. (2013) found that this formulation predicts satisfactorily the occurrence of blowing snow events at an alpine site over a 10 yr period provided the mechanical fragmentation of snow grains during blowing snow events is taken into account. U 5t is then converted into a threshold friction velocity u * th using the same method as Durand et al. (2005). Blowing snow occurs at grid points where the friction velocity u * is higher than u * th . u * is computed in SURFEX as a function wind speed at level of atmospheric forcing and transfer coefficient for momentum, C D : u * = √ C D U Surf . The expression of C D accounts for thermal buoyancy effects (Louis, 1979) and does not include particle buoyancy effects (Bintanja, 2000;Gallée et al., 2001).

Saltation layer
The saltation layer develops where snow transport occurs. In our model, the saltation layer contributes to the total snow transport and acts as a lower boundary condition for the suspension layer.
Sørensen (2004) proposed a physically-based formulation for the horizontal transport rate, Q Salt (kg m −1 s −1 ), of any particle in saltation: where V = u * /u * t . Constants a, b and c must be determined from observations (M. Sørensen, personal communication, 2012). For snow particles, we used the measurements of Nishimura and Hunt (2000) and found that a = 2.6, b = 2.5 and c = 2 allow to reproduce their observations. Doorschot and Lehning (2002) have similarly found that an earlier version of Eq. (11) (Sørensen, 1991) gives mass fluxes in a good agreement with those simulated by their saltation model and those observed by Nishimura and Hunt (2000). Q Salt is a stationary transport rate. Nemoto and Nishimura (2004) suggested that a time of 1-2 s is necessary to reach a steady state in the saltation layer. This corresponds to typical lengths of 1-20 m for wind speed close to the surface ranging from 1 to 10 m s −1 . These lengths are lower than the targeted horizontal resolution of Meso-NH/Crocus (50 m in this study, see Sect. 6.1) so that we can use a stationary mass flux. The model of Lehning et al. (2008) uses a similar assumption and manages to reproduce snow distribution down to an horizontal resolution of 5 m .
Mass exchange between the saltation layer and the suspension layer requires the computation of a reference concentration, c Salt (kg m −3 ), at the top of the saltation layer. The thickness of this layer is given by: h Salt = 0.08436u 1.27 * (Pomeroy and Male, 1992).  (Nishimura and Hunt, 2000): where λ is a dimensionless parameter (0.45 for snow, Nishimura and Hunt, 2000). u part denotes the horizontal particle velocity within the saltation layer which only depends on state of surface snow: u part = 2.8u * t (Pomeroy and Gray, 1990). The corresponding number concentration, N salt (m −3 ), is then computed from c salt assuming a two-parameter gamma distribution in the saltation layer and a fixed value for the average radius of particles in saltation, r mSalt . For the model application in alpine terrain, we determined r mSalt using data collected at Col du Lac Blanc by Naaim-Bouvet et al. (2011). The average radius in the saltation layer has not been measured but values in the suspension layer (between 25 and 29 cm) varies from 65 to 90 µm. Assuming that r m follows a power law r m = az b with b ≈ −0.25 (Schmidt, 1982), we found values for r mSalt ranging from 100 to 137 µm. We finally chose r mSalt = 110 µm. Note that the limitations of the formulations of the saltation layer are discussed in Sect. 7.2.

Snow erosion and accumulation
At each grid point, Meso-NH/Crocus computes a net mass flux, F Net (kg m −2 s −1 ), between the snowpack and the atmosphere: where F Salt and F Susp denote the contribution of the transport in saltation and in turbulent suspension, respectively. F Precip represents snowfall simulated by the micro-physical scheme of Meso-NH (Pinty and Jabouille, 1998). Figure 2 summarized the mass exchange between the snowpack and the atmosphere. The model of Lehning et al. (2008) simulates snow erosion and deposition in a similar way. This method lacks the feedback of the suspension on the saltation concentration. A prognostic equation for snow concentration in saltation would be required to overcome this limitation as in Gauer (1999). The current version of the model uses the steady state assumption for the saltation layer. The contribution from saltation corresponds to the divergence of the vector transport in saltation Q: where u MNH is the wind vector at the first level of Meso-NH. F Susp is a net mass flux between the saltation layer and the lowest level of the atmosphere. It follows: where F Sed is the sedimentation flux from the atmosphere and F Turb the turbulent flux of blown snow particles towards the atmosphere. F Sed corresponds to a loss of mass for the atmosphere: where c Surf denotes the near-surface concentration of blown snow particles (see Sect. 3.3) and V q the mass-weighted mean fall speed at this level. The second term of Eq. (15) F Turb follows the expression of Gallée et al. (2001): where U Surf denotes the near-surface wind speed (see Sect. 3.3). Similar fluxes are computed for the number concentration using N Salt and the near-surface number concentration N Surf (see Sect. 3.3). Erosion occurs where F Net < 0 and snow layers are removed from the snowpack profile simulated by Crocus. Therefore, snow layers with different characteristics may be exposed successively at the top of the snowpack during a blowing snow event. Snow accumulation is simulated where F Net > 0. Deposited snow is added to the existing snowpack using the routines of Crocus handling the layering of the snowpack. The current version of the model uses fixed values for the characteristics of deposited snow (density: 250 kg m −3 , dendricity: 0.3 and sphericity: 0.85). A future version of the model will account for the evolution of snow grain characteristics under snow transport.

Surface/atmosphere coupling
Mass fluxes between the snowpack and the atmosphere (Eqs. 16 and 17) require the estimation of near-surface variables (wind speed, U Surf , and number and mass concentration, N Surf and c Surf ). However, the vertical profile of blowing  Budd, 1966;Mann, 1998;Nishimura and Nemoto, 2005). 1-D blowing snow models (e.g. Déry and Yau, 1999;Mann, 1998) use a vertical grid with a very high resolution close to the ground to reproduce this gradient of concentration. For example, the model of Mann (1998) has a stretched grid with a first level at 10 cm above the surface and 70 levels in the lowest 100 m of the atmosphere. A similar configuration can be achieved with Meso-NH. However, such configuration is not suitable to simulate 3-D atmospheric flow in alpine terrain where the presence of large slopes generates large vertical velocities and requires small time steps to avoid numerical instabilities. An alternative solution is offered by the one-dimensional surface boundary layer (SBL) scheme Canopy (Masson and Seity, 2009). Canopy is implemented in SURFEX and includes prognostic atmospheric layers between the surface and the first level of Meso-NH. The evolution of SBL variables (wind, temperature, specific humidity and TKE) is resolved prognostically, taking into account large-scale forcing, turbulence, and, if any, drag and canopy forces and surface fluxes. These fluxes are computed between the surface and the lowest level of Canopy and sent back to the atmosphere. Therefore, Canopy allows to increase the vertical resolution near the ground in the surface scheme without the drawback of smaller time steps. As a consequence, blowing snow variables N s and q s have been implemented in Canopy to compute mass exchange between the snowpack and the atmosphere. The impact of this implementation is illustrated in Sect. 5.
The evolutions of N s and q s in Canopy are governed by the following one-dimensional equations: The terms in these equations are equivalent to the terms appearing in Eqs. (2) and (3). The contribution of advection at each Canopy level is computed through the advection in Meso-NH of the total mass and total number of blown snow particles in Canopy. This procedure is detailed in Appendix B. The sublimation terms S N and S q are computed as a function of temperature and specific humidity at each Canopy levels. Feedbacks on the SBL are included by adding the contribution of blowing snow sublimation in Canopy to the sensible and latent heat fluxes emitted towards the atmosphere. This modifies in return the temperature and humidity at the lowest levels of Meso-NH and, therefore, at Canopy levels.
Equations (18) and (19) are solved over the stretched grid of Canopy. For our application, it is made of 5 levels with a lowest level at 15 cm above the snowpack. When Canopy is activated, U Surf , N Surf and c Surf are taken at the lowest level of Canopy while the net mass flux at the top of Canopy, F Netκ , is sent to Meso-NH (Fig. 2).

Study site and evaluation strategy
Our study site is the experimental site of Col du Lac Blanc near the Alpe d'Huez ski resort, French Alps (Fig. 3). Due to the surrounding topography, the pass may be considered as a natural wind tunnel where 90 % of observed wind blows from two directions: northeast and south .
In the following sections, we propose a first evaluation of Meso-NH/Crocus using data collected at Col du Lac Blanc during two northern blowing snow events in 2011. Their main characteristics are given in Table 1. In-situ measurements collected during both events include: (i) meteorological conditions (wind speed and direction, air temperature) at three automatic weather stations (AWS) located around the pass, (ii) vertical profile (up to 3.5 m) of wind speed on a meteorological mast at the pass and (iii) vertical profiles of fluxes and radius of blown snow particles using three Snow Particles Counters (SPC, Sato et al., 1993) at the pass. Additionally, the evolution of snow depth was followed for the first event using data from a Terrestrial Laser Scanner (TLS, Prokop, 2008). Snow depths were measured over an area of 0.54 km 2 around the pass before (17 February 2011) and after (28 February 2011) the first event. Note that TLS measurements are not available for the second event.
The evaluation of Meso-NH/Crocus follows two steps. Results of 1-D simulations are firstly compared with vertical profiles of radius and fluxes of blown snow particles measured by the SPC during the 22-26 February event (Sect. 5). This allows us to discuss some of the model's features. 3-D simulations of the 18-19 March event are then presented and evaluated (Sect. 6). Figure 4 describes the meteorological conditions observed at Col du Lac Blanc for the 18-19 March event. 30 cm of fresh snow accumulated with a moderate wind during the night from 16 to 17 March have been redistributed by an intense wind blowing from the north. Light snowfall occurred on 19 March between 03:00 and 08:00. We selected the 18-19 March event as a case study for two main reasons: (i) the shorter duration of this event compared to the 22-26 February event (Table 1) makes it less computationally expensive to simulate in 3-D and (ii) the occurrence of snow transport without concurrent snowfall during most of the event allows us to focus the validation on the blowing snow scheme detailed in this paper. Simulation of blowing snow events with concurrent snowfall will be considered in future studies.  In the atmosphere the initial wind profile is logarithmic with a 5 m wind speed depending on the category and a roughness length z 0 = 1.3 × 10 −4 m. This value is based on observations collected at Col du Lac Blanc during northern blowing snow events (Vionnet, 2012). The initial nearsurface pressure is set to 730 hPa. The atmospheric boundary layer is stably stratified with potential temperature increasing gradually with height (vertical gradient: 0.005 K m −1 , Brunt-Väisälä frequency: 0.013 s) from a near surface value taken as the average potential temperature during the blowing snow event (288.5 K). Finally, the initial profile of relative humidity with respect to ice decreases logarithmically from 100 % at the surface to 50 % at 100 m. Above that level relative hu-midity is constant. At the bottom of the atmosphere the snow surface has a 5 m threshold wind speed equal to 10 m s −1 .
Three sets of 1-D simulations have been carried out with different vertical resolutions. Configuration LR (Low Resolution) uses a stretched vertical grid of 70 layers with 19 layers in the lowest 100 m of the atmosphere (lowest level: 3 m) while configuration HR (High Resolution) has 150 levels with 70 layers in the lowest 100 m of the atmosphere (lowest level: 0.15 m). Canopy is not activated in configurations LR and HR and the near-surface variables used in Eqs. (16) and (17) are taken at the first atmospheric level of Meso-NH. Finally, configuration LRC (Low Resolution Canopy) has the same vertical grid as configuration LR but uses Canopy which increases the vertical resolution close to the ground by adding 5 atmospheric layers (lowest Canopy level: 0.15 m). Each simulation last 20 min until stationary profiles are reached and blowing snow sublimation is not activated. The model time step is 0.1 s in configuration HR and 0.5 s in configuration LR and LRC.
Two additional sets of 1-D simulations in configuration LRC have been run to illustrate the impact of the parameterization of blowing snow sublimation. The first method (referred as DBL) computes the sublimation rate (Eq. 10) as a function of the PSD given by the double-moment scheme. The second set of simulations used the method described in GZ11 (referred as REP). A representative radius for blowing snow particles is used instead of the PSD simulated by the double-moment scheme. This representative radius  (Durand et al., 1993). corresponds to the single particle size that would give the same total sublimation rate for an equal concentration as an ensemble of particles based on a gamma size distribution with mean radius of 35 µm and α = 3. This results in a particle size of 62.5 µm. Note that this value is the same as the one used by GZ11 despite different values for the ventilation velocity and for α. Simulations including sublimation last 60 min. In configuration HR the model also reproduces the observed decrease in mean radius with height (Fig. 5, Left). This is a direct benefit from the double-moment scheme which allows size-sorting. However, simulated mean radius above 1 m are smaller than observations. This may be associated with the overestimation of low mean radius by the SPC which does not detect snow particles of radius lower than 25 µm (Nishimura and Nemoto, 2005). Note also that the model does not reproduce all the variability of mean radius close to the surface. This results from the assumption of a constant mean radius in the saltation layer.

Model results
The use of Canopy (configuration LRC) allows the model to simulate profiles of flux and radius similar to configuration HR close to the surface (0.15-2 m). Differences at higher levels are associated with the coarser vertical resolution in the range 2-10 m in configuration LRC. However, these differences have little consequences in terms of snow transport rate integrated over the lowest 100 m of the atmosphere (values Q HR and Q LRC reported on Fig. 5, Right). Therefore Canopy appears as an efficient tool to reproduce the strong gradient of blowing snow fluxes observed close to the surface. prognostically as in configuration HR but without the drawback of having a very fine vertical grid in the atmospheric model. The computational cost is divided by 8 in configuration LRC compare to configuration HR. Examples of local sublimation rates are shown on Fig. 6. They are maximum close to the ground for both methods (REP and DBL) which is a consequence of the large number of snow particles in this region. Method REP predicts higher sublimation rates than method DBL below 1 m. This results from the fact that method REP computes a sublimation rate representative of an ensemble of particle based on a gamma size distribution with r m = 35 µm while r m is higher than 35 µm below 1 m (Fig. 5, Left). Method DBL accounts for this effect on the sublimation rate while method REP does not. This has consequences on the hourly sublimation rate (vertically integrated over the lower 100 m of the atmosphere). It is multiplied by 2.20, 1.65 and 1.46 for 5 m wind speed of 12, 14 and 16 m s −1 , respectively, when using method REP instead of method DBL. Therefore, using a double-moment scheme allows to represent the observed decrease of mean radius with increasing height and to account for its consequences on the sublimation rate of blowing snow.

Model configuration
We carried out two 3-dimensional simulations for our case study: one with the sublimation of blowing snow (simulation SUBL) and one without (simulation CTRL). Both simulations started on 18 March at 01:00 and lasted until 18 March at 23:00. They covered the period of intense snow transport without concurrent snowfall (Fig. 4). The simulation domain was centred at Col du Lac Blanc and covered 3 × 3 km at an horizontal resolution of 50 m (Fig. 3). Meso-NH uses a vertical grid similar to the one used in configuration LRC (Sect. 5.1) with 20 layers in the lowest 200 m of the atmosphere. The height of the lowest level ranged from 1.9 to 3.1 m due to the terrain-following coordinates. Canopy added 5 extra atmospheric layers (lowest level: 15 cm above the snowpack). The snowpack was discretized in Crocus by 20 layers. The value for surface roughness was the same as the one used for 1-D simulations (z 0 = 1.3 × 10 −4 m). The model time step was 1.5 s.
The Meso-NH model was initialized on 18 March at 01:00 using a vertical atmospheric sounding. This sounding is first vertically interpolated on the model vertical grid and then distributed homogeneously over the simulation domain. A final adjustment allows the 3-D wind to fulfill the anelastic constraint and the free-slip boundary condition. This final wind field is used as initial field for the Meso-NH model. The vertical sounding is taken from a Meso-NH simulation at an horizontal resolution of 450 m at the grid point corresponding to Col du Lac Blanc. Input and forcing for this simulation comes from the analysis of the operational model AROME.
Near-surface wind in the sounding is adjusted to match the observed wind field at Col du Lac Blanc. This operation is repeated each hour to produce updated boundary conditions for the atmospheric model. This aims at reproducing the observed temporal evolution of wind speed at Col du Lac Blanc. This initialization method is only suitable for blowing snow event without snowfall. Accurately simulating cloud dynamics and resulting precipitations requires indeed grid-nesting techniques. The Crocus snowpack model was initialized using an horizontally homogeneous snowpit. 30 cm of erodible snow (U 5th = 9 m s −1 ) cover a layer of non-erodible snow as observed at the experimental site before the event. Such method allows to reduce model uncertainties about the initial snowpack but does not reproduce snowpack variability as a function of altitude, slope and aspect.

Atmospheric flow
The most sensitive transport modelling parameter in complex terrain is the driving wind field (e.g. Gauer, 1999;Lehning et al., 2008;Schneiderbauer and Prokop, 2011). In this section, we present an evaluation of the spatial and temporal evolution of the near-surface atmospheric flow simulated by Meso-NH. Figure 7a shows the modelled near-surface wind field at 07:00 on 18 March 2011 when maximum wind speed has been measured at Col du Lac Blanc (Fig. 4). Figure 7a illustrates how topography exerts a strong control on the wind direction and speed. The atmospheric flow is locally channeled along a north-south axis at Col du Lac Blanc (in between AWS Lac Blanc and Muzelle). The presence of ridges also modifies the flow (area A1 and A2 on Fig. 7a). Along theses ridges, the wind direction is normal to the ridge axis and the wind speed increases locally. Crest speedup is especially important in area A1 where the nearsurface wind speed reaches 33 m s −1 . We believe that the model tends to over-estimate the wind speed in this region.  identified similar over-estimation of crest speed-up in high-resolution (5 m) simulations. They associated it with an insufficient geometrical resolution of sharp crests in the model. Finally, the model simulates the formation of recirculation zones (area R1, R2 and R3 on Fig. 7a) in some leeward regions. In these areas, the wind direction is variable and may be opposed to the main flow while the wind speed is low. Therefore, Meso-NH is able to reproduce the main characteristics of wind flow in complex terrain : wind channelling, crest speedup and formation of recirculation zone downwind of the ridge. Such small scale flow features are the driving mechanisms behind the inhomogeneous snow distribution resulting from blowing snow events.
Wind speed and direction have been measured at three AWS around Col du Lac Blanc (Fig. 3). Table 2 shows that the model compares relatively well to the observations for the 2 AWS located on both side of the pass (Lac Blanc and Muzelle). Therefore, the radio-sounding method to initialise the atmosphere and to provide lateral boundary conditions allows the model to reproduce the evolution of wind speed at the pass during the blowing snow event. The main flow direction is also satisfactorily reproduced even if the model does not capture the eastern component of the wind at the Lac Blanc AWS. Results are different at the Dome AWS where the model represents well the wind direction but overestimates the wind speed. At this station, the atmospheric flow is highly turbulent compared to the two other stations. The averaged gust factors (defined as the ratio of maximum wind speed to average wind speed per period of 15 min) measured during the period of intense snow transport (02:00 to 18:00) are 2.08, 1.31 and 1.35 at AWS Dome, Muzelle and Lac Blanc, respectively. Meso-NH does not represent well this highly turbulent flow and predicts instead a flow with high wind speed. This will influence the intensity of snow redistribution.
The vertical profiles of wind speed close to the surface simulated by the model have been compared to wind speed profiles measured on a 4 m meteorological mast located at the pass (Fig. 7b). For this comparison, measured 15 min mean vertical profiles of wind speed were averaged into four 2 m wind speed bins (1 m s −1 ) ranging between 13 and 16 m s −1 . The same treatment has been applied to simulated profiles. Figure 7b shows that the model reproduces satisfactorily the near-surface vertical profile of wind speed for each category. Canopy turned out to be an adapted tool to refine the vertical profile of wind speed close to the snowpack. This is of prime importance to estimate with accuracy blowing snow fluxes.

Blowing snow fluxes
Meso-NH/Crocus simulates snow transport as a function of the type of surface snow and the wind forcing described at the previous section. Therefore, it is highly variable in space and time. Blowing snow fluxes have been measured at two levels above the surface (1.19 and 3.27 m on average) by two SPCs located at the pass. Fluxes at 1.19 m are roughly one order of magnitude higher than fluxes at 3.27 m (Fig. 8a)  evolution of blowing snow fluxes at both levels. Simulated snow transport stops at 17:00 when the wind speed drops below the threshold wind speed (Fig. 8b).
More quantitative results can be obtained by considering the total snow mass transported at both levels during our case study (Table 3). Meso-NH/Crocus tends to underestimate this snow mass at both levels (−11.0 % at 1.19 m and −44.4 % at 3.27 m). A potential explanation for this underestimation is the lower mean wind speed simulated by the model compared to observations ( Fig. 8b and Table 3). To better understand the model's behaviour as a function of the wind speed, the averaged simulated and observed blowing snow fluxes have been computed at both levels for five 2 m wind speed bins (1 m s −1 ) ranging between 12 and 16 m s −1 . Results are presented on Table 4. Meso-NH/Crocus simulates higher fluxes than the observations at 1.19 m. The relative bias between model and observations decreases when the wind speed increases (from 4.6 at 12 m s −1 to 1.2 at 15-16 m s −1 ). At 3.27 m, the model simulates higher fluxes in the range 12-14 m s −1 (mean relative bias of 2.3 for these 3 categories) and lower fluxes in the range 15-16 m s −1 (mean relative bias of 0.9 for these 2 categories). This analysis shows that model errors in term of total transported snow mass result from two main sources: (i) the simulation of the wind evolution during the blowing snow event and (ii) the blowing snow scheme itself which defines the quantity snow transported as a function of the wind speed. For our case study the underestimation of the total snow mass transport at both levels (Table 3) results from an underestimation of the wind speed at the location of the SPC. Figure 9 shows how snow has been redistributed in Meso-NH/Crocus during this blowing snow event. We present here a map of difference of snow water equivalent (SWE) between the end and the beginning of this event. The model simulates the formation of successive areas of snow erosion and deposition. This map is consistent with what is qualitatively expected from the wind field structure (Fig. 7a). Snow erosion is simulated on windward areas and deposition occurs in regions of decreasing wind speed. Areas of erosion are generally larger than the corresponding areas of deposition.

Patterns of snow erosion and deposition
Maps of snow depth difference measured by TLS measurements have been used in recent studies to evaluate the quality of the snow redistribution simulated by models Schneiderbauer and Prokop, 2011). However, this map is not available for our case study as mentioned in Sect. 4. Therefore, instead of a direct comparison, we propose here a qualitative analysis in term of typical patterns of erosion and deposition around Col du Lac Blanc for northern blowing snow events. For this purpose, we used the TLS data available for the 22-26 February event. This event presents similar conditions to our case study in terms of mean wind speed and direction (Table 1). Schirmer et al. (2011) have shown that individual storms arriving from the same direction produce similar patterns of snow depth changes at the Wannengrat catchment (Swiss Alps). The comparison presented in this section cannot be considered as a formal evaluation of the model ability to simulate snow redistribution but aims rather at exploring what is possible with the current model resolution of 50 m. Figure 10a shows the map of snow depth difference for the event of 22-26 February 2011. This map illustrates the presence of fine scale structures of erosion and deposition.
The Cryosphere, 8, 395-415, 2014 www.the-cryosphere.net/8/395/2014/ Erosion is mostly observed on the northern side of the pass where areas of snow deposition are also present. Their formation is related to topographic features such as cross-loaded slope (arrow 1) and slope change (arrow 2) responsible for the formation of a local recirculation zone on the windward side of the pass. On the southern side, deposition is generally observed with maxima of snow deposition at a steep slope break (arrow 3) and at the bottom of west-facing cliffs (arrow 4). Averaging to an horizontal resolution of 50 m (Fig. 10b) reduces greatly the variability of snow depth difference but preserves the general structure of erosion and deposition. Erosion is observed on the northern side of the pass while deposition occurred mainly on the southern side. Figure 10c shows that the model simulated erosion on northern side of the pass and deposition on the southern side for our case study. The location of the main deposition area (arrows 5, 6 and 7) is satisfactorily captured but their extension tends to be over-estimated (arrows 5 and 6). The formation of these area of deposition is generated by slope breaks that locally modified the atmospheric flow. A detailed analysis shows that these local topographic features are not present in the 50 m model orography. Therefore, the model at 50 m grid spacing reproduced satisfactorily the large-scale variability of snow depth difference around the pass (ridge scale) but missed smaller scale erosion and deposition patterns.

Influence of blowing snow sublimation
We carried out a second simulation of our case study including blowing snow sublimation (simulation SUBL). A comparison with the reference simulation (simulation CTRL) allows us to discuss the influence of blowing snow sublimation.
As explained in Sect. 3.1.2, blowing snow sublimation acts as a sink of snow mass and modifies in return near-surface meteorological variables. Figure 11 illustrates the effects of blowing snow sublimation on near-surface variables simulated by Meso-NH/Crocus at Col du Lac Blanc. Accounting for sublimation reduces suspended snow concentration of 4 % on average in the first 3 m of the atmosphere. During this event, sublimation rate ranges from 0 to −2.7 mm SWE day −1 . It is maximal at 07:00 when blowing snow concentration reaches its highest value. The evolution of near-surface variables also controls the sublimation rate. For a given snow concentration in the atmosphere, the decrease of humidity at the beginning of the event (03:00-09:00) changes the sublimation rate. It increases from −1 mm SWE day −1 at 04:00 to −2.5 mm SWE day −1 at 08:20 for the same snow concentration (1.5 g m −3 ) and a relative humidity dropping from 83 to 58 %. Sublimation effects on the surface boundary layer are illustrated in Fig. 11c and d. It leads to an increase of relative humidity with respect to ice. It is larger close to the surface with +9.5 % at 0.15 m and +3.6 % at 3.2 m on 18 March at 07:00. Blowing snow sublimation also affects potential temperature with a maximum decrease of −0.66 K at 0.15 m. The vertical gradient of potential temperature becomes stronger and increases the stability of the SBL at Col du Lac Blanc. Turbulent fluxes between the snowpack and the atmosphere are eventually modified. The detailed mass balance of a deposition area (Table 5) shows that the increase of relative humidity reduces surface sublimation by 29.7 %. Such reduction has been observed by Wever et al. (2009) in wind tunnel experiments. Overall, Meso-NH/Crocus simulates total sublimation (surface + blowing snow) loss of 1.46 mm for simulation SUBL and 0.47 mm for simulation CTRL. Therefore, for this case study and for the selected deposition area, the total sublimation is multiplied by 3 when accounting for blowing snow sublimation which becomes the main source of transfer of water vapour to the atmosphere (78 % of total sublimation).
Finally, blowing snow sublimation reduces snow concentration in the atmosphere and changes in snow redistribution are expected. Table 5 indicates that snow accumulation is reduced by 5.7 % for a deposition area covering 0.  the southern side of Col du Lac Blanc. Figure 12 confirms these results and shows that blowing snow sublimation reduces snow accumulation in areas of deposition identified on Fig. 9. The averaged deposited snow mass within our model domain is reduced by 5.3 % during this event in simulation SUBL compared to simulation CTRL. However, Fig. 12 also shows points where SWE at the end of the event is higher in simulation SUBL than simulation CTRL. To gain understanding on this unexpected increase, we computed the mass balance of a region located on the northern side of Col du Lac Blanc (green box on Fig. 12). It corresponds to an area of snow erosion (Fig. 9). Mean SWE difference between simulations SUBL and CTRL is +0.38 mm. 87 % of this increase is explained by a reduction of erosion flux resulting from a decrease of wind speed associated with higher atmospheric stability. The remaining 13 % are explained by the decrease of surface sublimation in simulation SUBL. Therefore, slight modifications of near-surface variables explain the increase of SWE observed at some points at the end of simulation SUBL compared to simulation CTRL.

Blowing snow sublimation
The influence of blowing snow sublimation has been estimated in a simulation including all sublimation feedbacks between snow and atmospheric dynamics. Results of this simulation can be compared with those from previous studies. Simulated sublimation rate range from 0 to −2.7 mm SWE day −1 . They are in the range of values measured by Pomeroy and Essery (1999) (−1.2 to −1.8 mm SWE day −1 ) and Schmidt (1982) (−0.51 and −5.27 mm SWE day −1 ). Blowing snow sublimation results in an increase of relative humidity of +9.5 % close to the surface. This increase is smaller than the maximum increase of 15 % simulated by GZ11 for a 43 h blowing snow event occurring over the Wannengrat catchment (2.4 km 2 , Swiss Alps). This difference may be explained by the formulation used by GZ11 which tends to give higher sublimation rates than the formulation used in Meso-NH/Crocus (Sect. 5.2). Note also that drier air was found for the period of maximum sublimation in their case study increasing blowing snow sublimation. Nonetheless, our results and those obtained by GZ11 are similar in a sense that they do not show saturation of the lowest levels of the atmosphere when blowing snow sublimation is activated. The significant increase of relative humidity near the snow surface with blowing snow sublimation mentioned by Déry et al. (1998) is not simulated in 3-D models in alpine terrain. This may be explained by the advection effects included in 3-D models and which are missing in 1-D unless explicitly included like in Bintanja (2001).
In simulation SUBL, we found that blowing snow sublimation causes a reduction of deposited snow mass of 5.3 % over the calculation domain. This value is higher than the 2.3 % found by GZ11 over the Wannengrat catchment. This difference can be explained by the meteorological conditions and the model formulation as mentioned before. The simulations also differ in terms of horizontal resolution (50 m in our study compared to 10 m in GZ11). Bernhardt et al. (2010) showed that their model predicts significantly smaller total amount of blowing snow sublimation over a winter season when decreasing horizontal grid spacing from 200 m to 30 m. This is due to the fact that crests where high sublimations rates are simulated have a smaller spatial extent at 30 m. Over our simulation domain, a finer horizontal resolution would probably result in a reduction of sublimation rates at the crests for the same reasons as in Bernhardt et al. (2010). However, quantifying this effect falls beyond the scope of this study and will require additional simulations that will be carried out in a future study.
Other studies focused on the impact of blowing snow sublimation over a whole winter at different spatial scales. For the Berchtesgaden park in Germany (210 km 2 ), Strasser et al. (2008) have found that 4.1 % of snowfall is lost by blowing snow sublimation. When including gravitational snow transport, the seasonal average loss is lowered to 1.6 % of annual snowfall for the same area (Bernhardt et al., 2012). In the Canadian Rocky Mountains, MacDonald et al. (2010) estimated that sublimation losses reach 17 to 19 % along a mountain crest (length: 210 m). Based on our results for a single event (reduction in snow deposition by 5.3 %), we can expect lower seasonal sublimation rates at Col du Lac Blanc since blowing snow events occur only during a fraction of the total time (10 % on average at Col du Lac Blanc, Vionnet et al., 2013). Groot Zwaaftink et al. (2013) confirmed this statement for the Wannengrat alpine catchment. They found that only 0.1 % of snowfall is lost by blowing snow sublimation over a winter while snow deposition can be reduced by 2.3 % for a single blowing snow event (GZ11).

Limitations in the coupled model
The fully coupled simulation of snow redistribution with Meso-NH/Crocus requires to run the model at a sufficiently high resolution for a time period covering the total duration of a blowing snow event (19 h on average at Col du Lac Blanc, Vionnet et al., 2013). In this study, we managed to simulate a 22 h blowing snow event in a fully coupled way at 50 m grid spacing. Results show that this configuration allows the model to capture the patterns of snow erosion and deposition at the ridge scale (Fig. 9). However, the detailed structures measured by TLS are not reproduced (Fig. 10). This constitutes the main limitation of the current version of the model. Indeed, the 50 m model orography does not represent local topographic features which modify the atmospheric flow and generate snow erosion and deposition.  found similar results at the Gaudergrat ridge (Swiss Alps) where their model at an horizontal resolution of 50 m reproduced patterns of snow deposition at the ridge scale but missed smaller scale deposition patterns. These patterns were partially reproduced using horizontal resolutions of 10 and 5 m. Using similar horizontal resolutions,  and Schneiderbauer and Prokop (2011) concluded that simulated snow redistribution compared satisfactorily with TLS measurements.
Capturing the small-scale erosion and deposition patterns will require to run Meso-NH/Crocus at 10 m grid spacing. Meso-NH can be used at this resolution in alpine terrain thanks to its 3-D turbulence scheme but there are some constraints linked to its dynamics. Due to the anelastic constraint, the pressure is diagnosed by solving an elliptic equation. In presence of steep orography, the convergence of the pressure solver is more difficult to reach. Specifically, the major limitation of this system is in terms of slope discontinuity: the topography used in the model must avoid a clifftype behaviour. The second constraint is linked to the eulerian numerical schemes, that limit the time step due to the Courant number. Amory (2012) ran Meso-NH at 12 m grid spacing around Col du Lac Blanc with a time step of 0.025 s. His simulation lasted 45 min which is not long enough to simulate the total duration of a blowing snow event. Furthermore, he showed that the border of the computational domain must be carefully chosen to avoid slope discontinuities. Therefore, the key issue to simulate snow redistribution at 10 m grid spacing lies in the ability of Meso-NH of running at this resolution for a time period covering a blowing snow event. This will be made possible by the on-going development of more efficient numerical schemes in Meso-NH. Model limitations have been identified in the formulation of the saltation layer. As detailed in Table 4  to overestimate blowing snow fluxes at 1.19 m. This may result from an overestimation of the mass concentration in the saltation layer which acts as a lower boundary condition for the suspension layer. The adaptation to snow of the formulation of Sørensen (2004) presented in Sect. 3.2.2 is based on a limited number of measurements and may not cover the whole range of threshold velocity. The formulation for the height of the saltation layer is also questionable. Meso-NH/Crocus uses the empirical formulation of Pomeroy and Male (1992) which tends to give lower values compared to the formulation used by Lehning et al. (2008) in alpine terrain (4.4 cm against 8.2 cm for example for u * = 0.6 m s −1 ). Using this formulation would decrease the concentration in the saltation layer for a given wind speed. The model also assumes that snow particle velocity in the saltation layer does not depend on the wind speed following the field measurements of Pomeroy and Gray (1990). However, more recent observations collected in wind tunnel showed that the particle velocity increases with increasing friction velocity for sand grains (Creyssels et al., 2009). Tominaga et al. (2013) showed that the ratio of snow particle velocity to wind velocity was about 40 % in the saltation layer. Using winddependent velocities for particles in the saltation layer would result in higher velocities and a decrease of snow concentration for a given mass flux. Overall, the current version of Meso-NH uses a simple description of the saltation layer which must be eventually replaced by a physically-based saltation model (e.g. Kok and Renno, 2009 The computation of blowing snow sublimation in Meso-NH/Crocus has also several limitations. Firstly, the model assumes that the ventilation velocity is equal to the settling velocity of suspended particles and neglects ventilation effects associated to turbulence (Dover, 1993;Bintanja, 2000). Furthermore, the effect of solar radiation is not included in the computation of the sublimation rate of an ice sphere. Measurements in a wind-tunnel (Wever et al., 2009) show that solar radiation may increase sublimation by 50 %. Finally, our model assumes that suspended snow particles are spherical. Based on the results of Thorpe and Mason (1966), Wever et al. (2009) estimated that the efficiency of the mass and energy exchange increased a factor 2.5 to 5 for dendritic shapes compared to hexagonal plates. The highest sublimation rates are observed with crystals composed of very fine branches. We estimate this effect is only active at the onset of blowing snow events before fragmentation modifies particles shape towards more rounded particles (Clifton et al., 2006;Vionnet et al., 2013). All these limitations lead to a probable underestimation of blowing snow sublimation rates by the model.

Conclusions
In this paper, we introduced the coupled system Meso-NH/Crocus dedicated to the study of snowpack/atmosphere interactions and wind-induced snow redistribution during blowing snow events. This system includes a blowing snow scheme which makes the distinction between snow transport in saltation and turbulent suspension and includes blowing snow sublimation. In the atmosphere, blown snow particles are represented by a double-moment scheme to capture the spatial and temporal evolution of the particle size distribution. At the surface, the model computes mass flux in saltation as a function of snow-surfaces properties given by Crocus and near-surface meteorological conditions. The resulting snow erosion or deposition includes the contributions of snow transport in saltation and turbulent suspension, and snowfall simulated by Meso-NH.
Meso-NH/Crocus has been evaluated in alpine terrain against data collected around the experimental site of Col du Lac Blanc (2720 m a.s.l., French Alps). Results of 1-D simulations indicate that the vertical resolution close to the ground must be sufficient to reproduce strong gradients of blowing snow concentration and compute mass exchange between the snowpack and the atmosphere. This is achieved in our model using a SBL scheme at the interface between Crocus and Meso-NH. The double-moment scheme for the suspension layer allows the model to capture satisfactorily the vertical profile of mean particle radius and to account for its consequence on sublimation rates of blowing snow.
A 22 h blowing snow event without concurrent snowfall was then selected and simulated in 3-D. At an horizontal resolution of 50 m, Meso-NH is able to capture the main features of near-surface atmospheric flow in alpine terrain including crest speed-up and the formation of recirculation zones. The total snow mass transported at two levels above the snowpack is underestimated by the model. This results from an underestimation of the wind speed at the location where blowing snow fluxes are measured. Simulated areas of erosion and deposition are consistent with what is qualitatively expected from the wind field structure: snow erosion on windward areas and deposition in regions of decreasing wind speed. However, at 50 m grid spacing, the model captures only the patterns of snow erosion and deposition at the ridge scale and misses smaller scale patterns observed by a terrestrial laser scanner around Col du Lac Blanc. Their formation is governed by topographic features of scale lower than 50 m.
In contrast with earlier-developed models, Meso-NH/Crocus includes most important feed-back mechanisms between snow and atmospheric dynamics, such as temperature and humidity effects of blowing snow sublimation. When activated for this case study, blowing snow sublimation reduces the total amount of deposited snow by 5 %. It also modifies the surface boundary layer in regions exposed to drifting snow and downwind of these zones with an increase of relative humidity by 5 to 9 % and a decrease of potential temperature by 0.3 to 0.7 K. Advection effects prevent the formation of a near-surface saturated layer. During this blowing snow event, blowing snow sublimation is the main source of water vapour transfer to the atmosphere and represents 78 % of total sublimation (blowing snow + surface).
Future developments will include the simulation of blowing snow events at higher horizontal resolution (up to 10 m) around Col du Lac Blanc to allow direct comparison with measurement from a terrestrial laser scanner. This will be made possible by the on-going development of more efficient numerical schemes in Meso-NH. A second application concerns the simulation of blowing snow events with concurrent snowfall which represent 37 % of observed blowing snow events at Col du Lac Blanc . The explicit representation of snowfall and blowing snow in the model will allow us to study the relative contribution of ground snow transport and preferential deposition of snowfall  to the spatial variability of the snowpack. This will require the use of grid-nesting techniques to produce atmospheric forcing from the synoptic scale down to the local scale.
Meso-NH/Crocus is the first fully coupled snowpack/atmosphere model that can simulate wind-induced snow transport at high-resolution in alpine terrain. It can also be applied over other snow-covered terrain such as ice sheets (Gallée et al., 2001) or over glaciers (Sauter et al., 2013). Its applications can be also extended to other topics regarding snowpack/atmosphere interactions such as the dynamics of katabatic wind over glacier (e.g. Claremar et al., 2012) or the influence of a patchy snow cover on the dynamics of the SBL