the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Understanding snow saltation parameterizations: lessons from theory, experiments and numerical simulations

### Daniela Brito Melo

### Armin Sigmund

### Michael Lehning

Drifting and blowing snow are important features in polar and high mountain regions. They control the surface mass balance in windy conditions and influence sublimation of snow and ice surfaces. Despite their importance, model representations in weather and climate assessments have high uncertainties because the associated physical processes are complex and highly variable in space and time. This contribution investigates the saltation system, which is the lower boundary condition for drifting and blowing snow models. Using a combination of (previous) measurements and new physics-based modeling with large-eddy simulation (LES), we show that the prevailing parameterizations that describe the saltation system in atmospheric models are based on contradictory assumptions: while some scaling laws are typical of a saltation system dominated by aerodynamic entrainment, others represent a saltation system controlled by splash. We show that both regimes can exist, depending on the friction velocity. Contrary to sand saltation, aerodynamic entrainment of surface particles is not negligible. It is important at low wind speeds, leading to a saltation height and near-surface particle velocity which increase with the friction velocity. In a splash-dominated saltation regime at higher friction velocities, the saltation height and near-surface particle velocity become invariant with the friction velocity and closer to what is observed with sand. These findings are accompanied by a detailed description of the theoretical, experimental and numerical arguments behind snow saltation parameterizations. This work offers a comprehensive understanding of the snow saltation system and its scaling laws, useful for both modelers and experimentalists.

- Article
(2386 KB) - Full-text XML
- BibTeX
- EndNote

The eolian or wind-driven transport of snow is common in snow-covered regions, in particular in high mountain areas, the Arctic and the Antarctic. This phenomenon, generally named drifting or blowing snow, occurs when the wind velocity is sufficiently high to lift the snow particles from the surface. During snow transport events, the wind erodes the snowpack in the most exposed regions and leads to snow re-deposition in sheltered areas. Therefore, the formation of bedforms, such as dunes, ripples, and sastrugi (Filhol and Sturm, 2015), as well as cornices at the leeward side of mountain crests (Hancock et al., 2020; Yu et al., 2023), is often observed. Snow sublimation is enhanced because the snow particles aloft are more exposed to the dry air flow (Dyunin, 1967; Liston and Sturm, 2004; Dai and Huang, 2014). In addition, some regions of the snow cover become harder and denser due to the deposition of snow particles that were previously aloft (Sommer et al., 2018). The influence of the eolian transport of snow on the mass and energy balances of the snow cover is therefore evident: it leads to snow redistribution; enhances the mass loss due to sublimation; and modifies the near-surface flow, temperature and moisture fields (Bintanja, 1998; Déry et al., 1998; Mott et al., 2018).

The terms drifting and blowing snow are commonly used to distinguish between weak and strong snow transport events. During drifting snow events, it is generally assumed that the snow particles do not rise more than 2 m above the ground, while during blowing snow events the snow grains are expected to reach higher regions of the atmosphere (e.g., Li and Pomeroy, 1997; Lenaerts et al., 2012). Nevertheless, in order to simplify the nomenclature, the term drifting snow is preferred throughout the text to designate snow transport events. These events are very frequent on the Antarctic continent, specifically in the coastal areas, which are subjected to strong katabatic winds throughout the year (Lenaerts et al., 2012; Palm et al., 2017). Drifting snow sublimation is also at a maximum along the coast due to the higher temperatures and lower values of relative humidity in comparison to the Antarctic Plateau (Gallée, 1998; Palm et al., 2017; Gerber et al., 2023). On average, the Antarctic surface mass balance is controlled by snowfall, but snow sublimation is the main mass ablation process (Lenaerts et al., 2012; Palm et al., 2017; Gerber et al., 2023). In some regions of East Antarctica, like Neumayer Station II and Adélie Land, drifting snow sublimation is expected to remove 17 % and 20 % of the snowfall, respectively (van den Broeke et al., 2010; Lenaerts and van den Broeke, 2012). Similar figures are obtained in the Low Arctic, such as in Trail Valley Creek, Canada, where drifting snow sublimation is expected to remove 19.5 % of the snowfall (Pomeroy et al., 1997). In some regions of the Canadian Prairies, estimates based on a simplified model indicate that sublimation losses due to snow transport may amount to 44 % or 74 % of the annual snowfall over a 4 km fetch depending on the meteorological conditions (Pomeroy et al., 1993). In the Alpine region, drifting snow sublimation is highly variable in space and is, on average, of the same order of magnitude of surface sublimation (Groot Zwaaftink et al., 2011). In mountain areas worldwide, snow transport and redistribution influence snow height variability, which plays a significant role in the occurrence of avalanches (Lehning et al., 2000; Das et al., 2012) and hinders the assessment of snowfall from direct snow height observations (Liston and Sturm, 2002; Mott and Lehning, 2010).

The current estimates of snow transport and sublimation are based on several decades of experimental measurements that supported the development of theoretical models and empirical correlations for the mass flow rate of particles aloft and their respective sublimation (e.g., Budd et al., 1966; Kobayashi, 1972; Schmidt, 1982; Pomeroy and Gray, 1990). Nevertheless, the overall contribution of drifting snow transport and sublimation to the surface mass balance of snow-covered regions on Earth is still uncertain. For example, while the Regional Atmospheric Climate Model (RACMO) (van Wessem et al., 2018) estimates a total of 102 ± 5 Gt yr^{−1} of drifting snow sublimation over the Antarctic continent (except the Antarctic Peninsula), the amount of drifting snow sublimation derived from satellite observations (Palm et al., 2017) over the Antarctic continent (except the region south of the 82° S parallel) yielded 393 ± 196 Gt yr^{−1}. The disagreement between the model and measurements of drifting snow sublimation brings into question the accuracy of large-scale predictions and highlights the need to revise the drifting snow formulations implemented in atmospheric models.

The motion of drifting snow particles is generally conceptualized in three modes: creep, saltation and suspension (Bagnold, 1941). Creep corresponds to the rolling and sliding of snow particles or to small hops with a maximum height of approximately the particle diameter. Saltation defines the ballistic motion of snow particles close to the surface, with a maximum height of no more than 10 to 15 cm (Gordon et al., 2009). These particles are either lifted by the turbulent flow (aerodynamic entrainment) or ejected by an impacting snow particle (splash). After impacting the surface, the saltating particles might deposit or bounce off into a new trajectory (rebound). This region close to the surface, where most particles are transported in saltation, is named the saltation layer. Suspension is the motion of small snow particles that travel without regular contact with the surface. These particles are advected by the wind for long distances and up to high regions of the boundary layer.

The general approach in atmospheric and snow models (Lehning et al., 2008; Vionnet et al., 2014; van Wessem et al., 2018; Amory et al., 2021; Sharma et al., 2023) is to parameterize the snow particle concentration and mass flux in the saltation layer and to compute the spatial and temporal distribution of snow particles in suspension by taking into account horizontal particle advection, particle settling and the turbulent vertical mixing of particles aloft (Dyunin and Kotlyakov, 1980; Pomeroy and Male, 1992; Déry et al., 1998). The bottom boundary of the snow suspension scheme is defined at the top of the saltation layer, where particle concentration is specified as a function of the near-surface wind velocity and some snow surface characteristics. The contribution of creep to the mass flux of snow particles is either neglected or assumed to be included in the snow saltation parameterization. Snow sublimation is usually estimated with the model of Thorpe and Mason (1966) (Groot Zwaaftink et al., 2011; Lenaerts et al., 2012; Vionnet et al., 2014; Sharma et al., 2023). The process of snow sublimation extracts heat and adds moisture from/to the near-surface flow, which is taken into account in the energy and moisture budgets of the near-surface atmosphere.

Shortcomings in the modeled drifting snow mass flux and sublimation can be related to uncertainties in the input parameters, such as the near-surface wind velocity, which greatly depends on the description of the surface topography (Gauer, 1998; Raderschall et al., 2008; Bernhardt et al., 2009; Lenaerts et al., 2012), or to inaccuracies in the drifting snow scheme itself. In fact, and concerning the latter, it seems that the resultant quantities of interest are significantly sensitive to the snow saltation parameterizations applied. For instance, using the Modèle Atmosphérique Régional (MAR), Amory et al. (2021) obtained an improved agreement between the simulated and measured drifting snow mass fluxes at D47 (Adélie Land, Antarctica) mainly by adjusting the minimum shear stress exerted by the flow that leads to the onset of snow drift. Another example can be given, this time with the atmospheric model RACMO: van Wessem et al. (2018) reduced the computed drifting snow sublimation over the Antarctic continent (except the Antarctic Peninsula) from 181 ± 9 to 102 ± 5 Gt yr^{−1} by halving the snow particle concentration at the top of the saltation layer. While the first parameter has the ability to vary the frequency of drifting snow events, the particle concentration at the top of the saltation layer mentioned in the second example affects the total mass of particles in suspension and, therefore, the rate of drifting snow sublimation.

As discussed above, a rigorous description of snow saltation in atmospheric models is of the utmost importance. Unfortunately, too little attention has been given to its development. The early model of Pomeroy and Gray (1990) has been used in the drifting snow schemes of RACMO (Lenaerts et al., 2012; van Wessem et al., 2018) and MAR (Amory et al., 2015, 2021). However, this saltation model is based on poorly constrained parameters and on the highly simplified assumption that the vertically averaged mass flux in saltation is well approximated by a point measurement inside the saltation layer. Vionnet et al. (2014) considered a different set of parameterizations in the atmospheric model Meso-NH. In particular, the vertical profile of particle mass flux in saltation is approximated by an exponential equation, as a function of the transport rate and the decay height, which we use to denote the decay parameter. This formulation was also adopted by Sharma et al. (2023) when coupling the Weather Research and Forecasting (WRF) model (Skamarock et al., 2019) with a drifting snow scheme and the snow model SNOWPACK (Bartelt and Lehning, 2002). The exponential decay of the mass flux in saltation with increasing height above the surface is a well-known feature of the saltation system (e.g., Sugiura et al., 1998; Nishimura and Nemoto, 2005; Kok and Renno, 2009; Martin and Kok, 2017). However, the calculation of the transport rate and the decay height is still a challenge. Several expressions have been proposed, but there is still no full consensus on how these quantities vary with the wind field and the snow cover characteristics.

In this work, we revisit the different snow saltation parameterizations available in the literature and their use in atmospheric models. In particular, we focus on the theoretical, experimental and numerical arguments that support them (Sect. 2). We highlight where consensus and disagreement exist in order to help the reader develop a critical view on the currently used snow saltation parameterizations. By pinpointing where most of the uncertainty lies, this exposition also seeks to motivate further studies in the field. In addition, we use a numerical model to investigate some aspects of the saltation system that are not yet fully understood (Sect. 3). The numerical model computes the particle trajectories in a turbulent flow by means of a large-eddy simulation flow solver coupled with a Lagrangian stochastic model (LES-LSM) (Sharma et al., 2018; Comola et al., 2019b; Sigmund et al., 2022; Melo et al., 2022). The transport rate and the vertical profiles of particle mass flux and velocity are investigated, as well as the particle hop height and length during saltation. Special focus is given to the decay height and to its relation to the saltation dynamics. In Sect. 4, we summarize the main results and give suggestions for future work.

In this section, the different snow saltation parameterizations available in the literature and their use in atmospheric models are presented and discussed. We focus on the most common approximations within the steady-state framework and highlight the evidence that support them. After an introduction about the steady-state assumption (Sect. 2.1), the parameterizations used to compute the transport rate (Sect. 2.2) and the vertical profiles of particle concentration, streamwise velocity and mass flux (Sect. 2.3) are analyzed. These quantities are needed to compute the particle concentration at the top of the saltation layer – the lower boundary condition for the advection and diffusion of suspended particles (Lehning et al., 2008; Vionnet et al., 2014; van Wessem et al., 2018; Amory et al., 2021; Sharma et al., 2023). In particular, the effect of wind velocity on the different quantities of interest is discussed. The saltation dynamics is also highly dependent on the granular bed characteristics, which are of particular importance in the study of snow saltation (Comola and Lehning, 2017). However, considerably less literature is available on this topic.

## 2.1 The steady-state assumption

Let us denote by *x* [m], *y* [m] and *z* [m] the streamwise, crosswise and vertical directions, respectively, in a Cartesian coordinate system. In this reference frame, we consider a turbulent air flow driven by a constant streamwise pressure gradient over a flat, homogeneous and erodible surface. It is assumed that the wind field is statistically steady and that the time-averaged wind velocity does not vary along the streamwise and crosswise directions. Under these conditions, the saltation system is considered statistically steady and invariant along the horizontal directions (steady-state saltation).

Close to the surface and in the absence of saltating particles, the time-averaged fluid shear stress can be assumed constant in height (Prandtl, 1935). This region is generally named the inner layer. In the atmospheric boundary layer, this constant shear stress assumption is valid in approximately the first 10 % of the boundary layer depth. In geophysical flows, this region is called the surface layer and has a thickness of 10–100 m (Stull, 1988). When the surface is aerodynamically smooth, this near-surface region is conceptualized in three sublayers: a viscous sublayer in the vicinity of the surface, a logarithmic sublayer in the upper region and a buffer sublayer in between. When the surface is fully rough (which is frequently the case in atmospheric boundary layer flows), the viscous sublayer is disrupted by the roughness features (Monin, 1970). In this case, the constant shear stress region is better described in two sublayers: the roughness sublayer and the logarithmic sublayer (Jiménez, 2004). In the latter, the turbulent shear stresses (often denoted as Reynolds stresses) are dominant in comparison to those of viscous origin. In addition, under neutral stability conditions, the average streamwise wind velocity is expected to follow a logarithmic profile, which is a function of the friction velocity, *u*_{∗} [m s^{−1}], and of the roughness length, *z*_{o} [m], which characterizes the roughness of the surface. In this way, the time-averaged fluid shear stress in the logarithmic sublayer, $\mathit{\tau}={\mathit{\rho}}_{\mathrm{a}}{u}_{\ast}^{\mathrm{2}}$ [Pa], where *ρ*_{a} [kg m^{−3}] is the air density, can be computed from the covariance of the velocity fluctuations or estimated by fitting the vertical profile of average streamwise wind velocity to a logarithmic function. In atmospheric models, the streamwise wind velocity is assumed to follow a logarithmic function between the surface and the first grid point (Gallée et al., 2001; White, 2003; Skamarock et al., 2019). In agreement with the Monin–Obukhov similarity theory (Monin and Obukhov, 1954), stability corrections are added to the standard logarithmic profile, which changes slightly the relationship between *u*_{∗}, *z*_{o} and the streamwise wind velocity at the first grid point above the surface. Under the constant shear stress assumption, the time-averaged fluid shear stress in the logarithmic sublayer is equal to the time-averaged surface shear stress, ${\mathit{\tau}}_{\mathrm{s}}={\mathit{\rho}}_{\mathrm{a}}{u}_{\ast ,\mathrm{s}}^{\mathrm{2}}$ [Pa] (henceforth referred to as surface shear stress), where ${u}_{\ast ,\mathrm{s}}$ [m s^{−1}] is the surface friction velocity.

When surface particles are entrained by the fluid flow and advected by it, there is an exchange of momentum between the fluid and the particles aloft. This leads to a decrease in the fluid shear stress in the saltation layer (Owen, 1964; Raupach, 1991). Above the saltation layer and within the logarithmic sublayer, the time-averaged fluid shear stress is still well approximated by *τ* and the average streamwise wind velocity profile is still expected to follow a logarithmic trend. However, this logarithmic profile is no longer a function of *z*_{o}, but of an equivalent roughness length, greater than *z*_{o}, which takes into account the momentum deficit in the saltation layer (Bagnold, 1941; Owen, 1964). At the surface, the fluid shear stress, *τ*_{s}, reaches a value smaller than *τ*.

During steady-state saltation, the net erosion and deposition of particles from/on the granular bed is equal to zero. In this way, for a given friction velocity, *u*_{∗} (characteristic of the logarithmic region outside the saltation layer), the rate at which particles are transported by the wind stays at an equilibrium value. The saltation system is therefore characterized by a type of self-regulatory behavior: a local increase in the particle flow rate above its equilibrium value locally extracts more momentum from the wind, reducing its velocity; this decreases the ability of the air flow to accelerate the particles in saltation and, consequently, decreases the particle flow rate back to its equilibrium value (e.g., Bagnold, 1941; Owen, 1964).

The premise behind the steady-state saltation assumption contrasts with the reality of snow-covered regions, which are rarely flat or homogeneous. In addition, the wind field in a mountainous terrain can seldom be characterized as statistically steady. Nevertheless, steady-state saltation is frequently assumed, based on the notion that the saltation system can adapt very rapidly to changes in the wind velocity (e.g., Anderson et al., 1991). Wind tunnel and field measurements have revealed a strong coupling between the saltation system and the wind velocity fluctuations (Paterna et al., 2016; Aksamit and Pomeroy, 2018). However, atmospheric models require simple and computationally light models that cannot capture the full complexity of the wind–particle interaction. In addition, the timescales of atmospheric models do not encompass the full wind field variability. In this context, when deriving snow saltation parameterizations for large-scale models, the assumption of steady-state saltation is regarded as a reasonable approximation.

## 2.2 The transport rate

The transport rate of saltating particles, *Q* [$\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$], corresponds to the mass flow rate of particles moving in saltation per unit width. From its definition, *Q* can be computed by multiplying the mass of particles in saltation per unit surface area by the average particle streamwise velocity.

The transport rate can also be computed as a function of the wind field characteristics by performing a balance of linear momentum in the streamwise direction to the saltating particles and to the fluid flow. Let us assume that the only force applied to a given saltating particle in the streamwise direction is the aerodynamic force (neglecting, for instance, electrostatic forces and mid-air collisions). In this case, the average aerodynamic force in the streamwise direction applied to a particle along a hop equals the average rate of change in the particle streamwise momentum, $m({v}_{x}^{\mathrm{i}}-{v}_{x}^{\mathrm{e}})/{t}_{\mathrm{p}}$, where *m* [kg] is the particle mass; ${v}_{x}^{\mathrm{i}}$ [m s^{−1}] and ${v}_{x}^{\mathrm{e}}$ [m s^{−1}] are the streamwise components of the particle velocity at impact and ejection from the surface, respectively; and *t*_{p} [s] is the duration of the particle hop (Fig. 1a). The balance of linear momentum in the streamwise direction applied to the fluid flow is simplified by assuming steady-state saltation conditions and by neglecting the pressure gradient. According to Newton's third law, the aerodynamic force applied to a particle is equal in magnitude and opposite to the force applied to the fluid. Therefore, the sum of aerodynamic forces in the streamwise direction applied to all particles in saltation per unit surface area equals the difference in shear stresses, *τ*−*τ*_{s} (Fig. 1b). From these considerations and assuming that the variety of particle hops can be approximated by an average trajectory (characteristic path) (e.g., Bagnold, 1941; Ungar and Haff, 1987; Durán et al., 2011; Kok et al., 2012), the following expression for the transport rate is obtained:

where *v*_{p} [m s^{−1}] is the hop-averaged particle streamwise velocity, $\mathrm{\Delta}{v}_{x,\mathrm{s}}={v}_{x}^{\mathrm{i}}-{v}_{x}^{\mathrm{e}}$ [m s^{−1}] is the difference between the impact and ejection streamwise velocities, and *l*_{p}=*v*_{p}*t*_{p} [m] is the particle hop length in the streamwise direction (henceforth, hop length). The angle brackets, 〈〉, represent averages over the ensemble of particles in saltation. Therefore, 〈*v*_{p}〉 denotes the average particle streamwise velocity. When deducing Eq. (1), the snow surface was assumed homogeneous and erodible. However, the presence of non-erodible roughness elements can be taken into account by adding a parameterization for shear stress partition (Pomeroy and Gray, 1990; Raupach et al., 1993).

The use of average quantities to describe the various particle trajectories introduces an error in the calculation of the transport rate. This assumption is not sufficiently discussed in the literature (e.g., Shao and Mikami, 2005). Nevertheless, Eq. (1) is the basis for several simple parameterizations (e.g., Bagnold, 1941; Ungar and Haff, 1987; Pomeroy and Gray, 1990; Sørensen, 1991; Kok et al., 2012) and is expected to yield a reasonable approximation for the scaling of *Q* with *u*_{∗}. The variety of parameterizations proposed for *Q* arise from different expressions for *τ*_{s}, 〈Δ*v*_{x,s}〉 and 〈*l*_{p}〉. These quantities are discussed in the next subsections.

### 2.2.1 The surface shear stress

The onset of saltation occurs when the surface shear stress grows above a given threshold, which allows the aerodynamic entrainment of surface particles to take place (Bagnold, 1941). This threshold shear stress is named the fluid threshold, *τ*_{ft} [Pa]. By definition, ${\mathit{\tau}}_{\text{ft}}={\mathit{\rho}}_{\mathrm{a}}{u}_{\ast ,\text{ft}}^{\mathrm{2}}$, where ${u}_{\ast ,\text{ft}}$ [m s^{−1}] is the fluid threshold friction velocity. According to Owen (1964), during steady-state saltation, the particles are continuously set in motion by aerodynamic entrainment, and the particle bombardment at the surface is only responsible for dislodging or loosening the grains. Naturally, he claims that saltation can only be sustained if *τ*_{s} is equal to or greater than a minimum value that allows the entrainment of loose grains, previously dislodged by impacting particles. Owen (1964) goes one step further and hypothesizes that the surface shear stress stays at this minimum value at which saltation is sustained (Owen's hypothesis): if *τ*_{s} is lower than this threshold, no aerodynamic entrainment of previously dislodged particles occurs, and, if it is larger than that threshold, more particles are entrained, which ultimately restores *τ*_{s} to its original value. According to Owen, the minimum surface shear stress required to sustain saltation is equal to the impact threshold, ${\mathit{\tau}}_{\text{it}}={\mathit{\rho}}_{\mathrm{a}}{u}_{\ast ,\text{it}}^{\mathrm{2}}$ [Pa]: the lower fluid shear stress outside the saltation layer for which saltating particles maintain a steady-state motion (Bagnold, 1941). ${u}_{\ast ,\text{it}}$ [m s^{−1}] is the impact threshold friction velocity. Due to the contribution of the impacting particles to dislodging the grains, Owen (1964) expects the impact threshold to be smaller than the fluid threshold. Indeed, Bagnold (1941) showed with wind tunnel experiments developed with sand that the shear stress needed to keep saltating particles in motion is smaller than that required to initiate saltation from rest. In his experiments, the shear stress was computed from fitting the wind velocity outside the saltation layer to a logarithmic function. Therefore, the reported values of shear stress are characteristic of the logarithmic sublayer. Nevertheless, at the fluid and impact thresholds, the surface shear stress and the shear stress outside the saltation layer should assume similar figures due to the low number of particles aloft and the consequent small momentum exchange between the fluid and the particles.

The experimental assessment of Owen's hypothesis is challenging because it requires high-frequency measurements of the wind velocity close to the surface within a saltation layer (e.g., Li and McKenna Neuman, 2012). The measurements of Walter et al. (2014) of the surface shear stress during saltation revealed a slight variation in *τ*_{s} as *u*_{∗} increases. In addition, some numerical models that include a statistical description of rebound and splash obtained a significant variation in the surface shear stress with respect to *u*_{∗} (Werner, 1990; Nemoto and Nishimura, 2004; Kok and Renno, 2009), while others predicted a negligible change (Niiya and Nishimura, 2017; Melo et al., 2022). Doorschot and Lehning (2002), who modeled steady-state saltation with a simplified numerical model, also obtained a slight variation in *τ*_{s} with the rise in *u*_{∗}. From the different trends obtained by measurements and models, it is clear that the self-regulatory behavior of steady-state saltation is not necessarily characterized by a constant surface shear stress, invariant with respect to *u*_{∗}, but that different equilibrium states can emerge depending on the wind and particle characteristics. In this way, Owen's hypothesis is expected to oversimplify the saltation dynamics (Kok et al., 2012). Nonetheless, the assumption that *τ*_{s} does not vary with *u*_{∗}, but depends solely on the bed characteristics, is regarded as an acceptable first-order approximation (Walter et al., 2014; Niiya and Nishimura, 2022). Therefore, it might be appropriate for simple models and parameterizations, regardless of aerodynamic entrainment being considered the main entrainment process (e.g., Pomeroy and Gray, 1990; Sørensen, 1991, 2004).

The numerical model of Kok and Renno (2009) has revealed an increase in *τ*_{s} with the increase in the particle diameter of the granular bed (see Fig. 13 in Kok et al., 2012). The same is found by Melo et al. (2022) when considering uniformly sized particles. In the latter, it is seen that other characteristics of the granular bed, as the particle size distribution and the existence of cohesive bonds between the particles, also influence the value of *τ*_{s} during steady-state saltation.

^{∗} Some authors define the impact threshold as the minimum surface shear stress for which saltation is sustained by splash (Comola et al., 2022). This definition of *τ*_{it} will be considered in Sect. 3.3.

In the drifting snow schemes of the atmospheric models RACMO, MAR, Meso-NH and CRYOWRF (Lenaerts et al., 2012; Amory et al., 2021; Vionnet et al., 2014; Sharma et al., 2023), *τ*_{s} is assumed invariant with respect to *u*_{∗} and is approximated by the fluid threshold, *τ*_{ft}. In these models, the fluid threshold varies with the snow surface characteristics, which change over time and space according to the meteorological conditions (Guyomarc'h and Mérindol, 1998; Gallée et al., 2001; Lehning et al., 2000). Therefore, the equality between *τ*_{ft} and *τ*_{s} guarantees that *τ*_{s} varies with the snow surface characteristics as proposed by Owen (1964). However, the assumption that *τ*_{s}=*τ*_{ft} neglects the effect of particle impacts in dislodging the grains, which justifies Owen's claim that *τ*_{s} must be smaller than *τ*_{ft}. The assumption that *τ*_{s}=*τ*_{ft} also implies that the fluid shear stresses for saltation initiation and cessation are equal, which is a useful assumption (note that Eq. (1) reduces to zero when *τ* equals *τ*_{s}). However, it is neither in agreement with wind tunnel measurements nor saltation models developed for sand (e.g., Bagnold, 1941; Kok et al., 2012). In the case of snow saltation, the relationship between *τ*_{s} and *τ*_{ft} is still not clear. While some numerical models predict *τ*_{s}>*τ*_{ft} (Nemoto and Nishimura, 2004; Melo et al., 2022), others predict *τ*_{s}<*τ*_{ft} (Niiya and Nishimura, 2017). Both trends are found with the model of Doorschot and Lehning (2002) depending on the characteristics of the snow cover. In this way, further investigation is needed to better assess the value of *τ*_{s} during snow saltation and the impact of the currently used assumption on the accuracy of the transport rate. A summary of the different shear stresses defined in this work and the main assumptions regarding the surface shear stress are presented in Table 1.

### 2.2.2 The ejection and impact velocities

When assessing the transport rate, Bagnold (1941) assumed that the ratio $\langle {l}_{\mathrm{p}}\rangle /\langle \mathrm{\Delta}{v}_{x,\mathrm{s}}\rangle $ scales with $\langle {v}_{z}^{\mathrm{e}}\rangle /g$, where ${v}_{z}^{\mathrm{e}}$ [m s^{−1}] is the vertical component of the ejection velocity and *g* [m s^{−2}] is the acceleration of gravity. One arrives at this result by solving the ballistic motion of a particle in saltation subjected to some simplifying assumptions (see, for instance, Sørensen, 1991). In addition, Bagnold (1941) assumed that $\langle {v}_{z}^{\mathrm{e}}\rangle $ is proportional to *u*_{∗}, based on the idea that the particle velocity at impact must scale with *u*_{∗} (please note that *u*_{∗} is the friction velocity outside the saltation layer and not the surface friction velocity, which is defined by ${u}_{\ast ,\mathrm{s}}$).

Numerical models, like the one proposed by Ungar and Haff (1987), generally suggest that 〈Δ*v*_{x,s}〉 is invariant with respect to *u*_{∗}. Ungar and Haff (1987) arrived at this result by assuming that in steady-state saltation each particle impacting the surface is replaced by a new one, equal in size, splashed from the surface upon impact. Indeed, during steady-state saltation, the zero net erosion and deposition of particles implies an equality between the mass flow rate of particles that fail to rebound and the mass flow rate of particles that successfully enter saltation via splash or aerodynamic entrainment (Kok et al., 2012). Contrary to Owen (1964), most numerical models assume that the surface processes of rebound and splash play an important role in steady-state saltation (e.g., Kok and Renno, 2009; Durán et al., 2011). These models compute the probability of rebound and the number and initial velocity of splashed grains as a function of the bed characteristics and the impacting particle mass and velocity. As a result, they predict that 〈Δ*v*_{x,s}〉 does not vary with *u*_{∗}. This assumption is present in the analytic model of Sørensen (2004), which is used in Meso-NH (Vionnet et al., 2014) and CRYOWRF (Sharma et al., 2023) to compute the transport rate. In addition, numerical models of snow saltation that represent the surface processes of aerodynamic entrainment, rebound and splash have also revealed ejection and impact velocities that are mainly invariant with respect to *u*_{∗} (Melo et al., 2022; Niiya and Nishimura, 2022).

Field and wind tunnel measurements conducted with sand show a negligible influence of *u*_{∗} on the near-surface particle velocity. For instance, Namikas (2003) found a better agreement between model results and field measurements of particle mass flux when assuming an averaged ejection velocity invariant with respect to *u*_{∗} for friction velocities ranging from 0.27 to 0.63 m s^{−1}. Moreover, wind tunnel experiments developed by Creyssels et al. (2009) and Ho et al. (2011) over an erodible sand surface show a negligible variation in the particle streamwise velocity close to the surface with respect to *u*_{∗}, considering friction velocities ranging from 0.24 to 0.67 m s^{−1} and from 0.42 to 1 m s^{−1}, respectively. During the experiments, the particles were illuminated by a laser sheet and recorded with a high-speed camera. The particle velocities were calculated with the particle tracking velocimetry (PTV) technique, and the particle streamwise velocity close to the surface was extrapolated from a linear regression of the average streamwise velocity profiles in the first 5 and 3 cm above the surface, respectively. The particle streamwise velocity close to the surface is named slip velocity and represents the average between the ejection and impact streamwise velocities.

The streamwise velocity of saltating snow was investigated by Aksamit and Pomeroy (2016) during field experiments in Alberta, Canada, using a similar experimental setup to the one used by Creyssels et al. (2009) and Ho et al. (2011). The streamwise velocity of ascending particles in the first 2.5 cm above the surface was reported for friction velocities ranging from 0.21 to 0.54 m s^{−1} for three different snow conditions (old snow, fresh snow and fine decomposing grains). The average particle ejection velocity in the streamwise direction, $\langle {v}_{x}^{\mathrm{e}}\rangle $, was extrapolated from a linear regression of the average streamwise velocity of ascending particles along the first 1 cm above the surface Even though steady-state conditions were not attained during the experiments, different trends were obtained for different snow surface characteristics. For the old and fresh snow covers, the streamwise ejection velocity was not found to vary monotonically with the friction velocity, while for the one composed of fine decomposing grains, a slight increase in *u*_{∗} was observed. Among the three snow surfaces, the latter was characterized by the highest hand hardness index. The increase in the streamwise ejection velocity with *u*_{∗}, obtained over the surface of fine decomposing grains, is in agreement with the increase in the particle slip velocity as *u*_{∗} increases, obtained by Ho et al. (2011) after performing wind tunnel measurements of saltating sand over a rigid (non-erodible) surface, considering friction velocities ranging from 0.3 to 0.45 m s^{−1}. Aksamit and Pomeroy (2016) obtained higher values of $\langle {v}_{x}^{\mathrm{e}}\rangle $ over the fresh snow cover than over the fine decomposing grains and old snow covers. In contrast, Ho et al. (2011) obtained higher values for the particle slip velocity over the rigid surface in comparison to the erodible surface.

In summary, most of the available literature suggests that the average ejection and impact velocities of particles in saltation are invariant with respect to *u*_{∗}. However, the few experiments developed with snow reveal different trends depending on the snow surface characteristics. In order to understand the particular case of snow saltation, more experimental measurements developed with snow are required, which must include a detailed characterization of the snow cover, the wind field and the kinematics of the particles in saltation.

### 2.2.3 The hop length

The particle hop length, *l*_{p}, is equal to the duration of the particle hop, *t*_{p}, times the hop-averaged particle streamwise velocity, *v*_{p}. If we approximate the variety of particle trajectories by an average trajectory, the average hop length, 〈*l*_{p}〉, can be estimated from the product 〈*t*_{p}〉〈*v*_{p}〉. 〈*t*_{p}〉 is expected to deviate from the theoretical value of $\mathrm{2}\langle {v}_{z}^{\mathrm{e}}\rangle /g$ (obtained when neglecting the aerodynamic forces in the particle motion), but it is nonetheless expected to scale with $\langle {v}_{z}^{\mathrm{e}}\rangle $
(Sørensen, 1991). In addition, as discussed in the previous section, the particle ejection velocity during steady-state saltation is generally regarded as invariant with respect to *u*_{∗}. Under these conditions, the scaling of 〈*l*_{p}〉 with *u*_{∗} is mainly dependent on the scaling of 〈*v*_{p}〉 with *u*_{∗}. Assuming that the particle ejection velocity does not vary with *u*_{∗}, the average particle streamwise velocity, 〈*v*_{p}〉, depends solely on the aerodynamic forces applied on the particles during their trajectories and, therefore, on the wind velocity inside the saltation layer.

Ungar and Haff (1987) concluded that 〈*l*_{p}〉 is invariant with respect to *u*_{∗} and proportional to the particle diameter, *d* [m]. This result is obtained by assuming that the ejection and impact velocities scale with (*d**g*)^{½} and that the wind velocity in the saltation layer is invariant with respect to *u*_{∗}. In contrast, the analytic model of Sørensen (2004), used in Meso-NH (Vionnet et al., 2014) and CRYOWRF (Sharma et al., 2023), arrived at an average particle hop length proportional to *u*_{∗}. In the model, the particle ejection velocity is also assumed invariant with respect to *u*_{∗}, and the logarithmic (particle-free) wind profile is modified to take into account the reduction in fluid shear stress induced by the particles in saltation. Therefore, the wind velocity decreases in comparison to the logarithmic profile, but it is nonetheless proportional to *u*_{∗} inside the saltation layer.

Wind tunnel measurements show that the wind velocity during saltation is expected to be invariant with respect to *u*_{∗} in the first 1 cm above the surface, where most saltating particles are (Bagnold, 1941; Li and McKenna Neuman, 2012). However, at higher elevations, the wind velocity is expected to scale with *u*_{∗} (e.g., Nishimura et al., 2014; Aksamit and Pomeroy, 2016). Therefore, if a significant number of particles travel above a height of 1 cm during a significant extent of their trajectories, 〈*v*_{p}〉 might increase as *u*_{∗} increases. This is in agreement with the field experiments of Namikas (2003) developed with sand, which revealed a slight increase in the average hop length with *u*_{∗}.

The numerical model of Durán et al. (2011) suggests that the scaling of the hop length with the friction velocity differs according to the saltation regime. At a low-velocity regime (for, approximately, ${u}_{\ast}/{u}_{\ast ,\text{it}}<\mathrm{4}$), most saltating particles are within this first 1 cm above the surface where the wind velocity is invariant with respect to *u*_{∗}. Therefore, the average particle velocity and hop length are expected to be invariant with respect to *u*_{∗}. However, at a high-velocity regime (${u}_{\ast}/{u}_{\ast ,\text{it}}>\mathrm{4}$), a considerable number of saltating particles reach higher elevations above the surface, where the wind velocity increases with *u*_{∗}. Under these conditions, both the average particle streamwise velocity and hop length are expected to increase linearly with *u*_{∗}. Taking into account that ${u}_{\ast}/{u}_{\ast ,\text{it}}>\mathrm{4}$ rarely occurs in nature in the case of saltating sand, Kok et al. (2012) argued that neglecting the effect of *u*_{∗} on the hop length can be regarded as a reasonable approximation. If both 〈Δ*v*_{x,s}〉 and 〈*l*_{p}〉 are assumed invariant with respect to *u*_{∗}, the transport rate scales with ${u}_{\ast}^{\mathrm{2}}$ (note that $\mathit{\tau}=\mathit{\rho}{u}_{\ast}^{\mathrm{2}}$) (Ungar and Haff, 1987; Durán et al., 2011). Conversely, if 〈*l*_{p}〉 increases linearly with *u*_{∗}, one arrives at $Q\propto {u}_{\ast}^{\mathrm{3}}$ (Sørensen, 1991, 2004).

Further studies are needed in order to fully understand the scaling of 〈*l*_{p}〉 with the friction velocity during snow saltation. In particular, it would be interesting to unveil if the two regimes proposed by Durán et al. (2011) are also representative of saltating snow. The full vertical profile of the particle streamwise velocity is discussed in Sect. 2.3.2.

### 2.2.4 The average particle streamwise acceleration

The hop-averaged particle streamwise acceleration, *a*_{p} [m s^{−2}], is given by the ratio $\mathrm{\Delta}{v}_{x,\mathrm{s}}/{t}_{\mathrm{p}}$. When assuming a characteristic particle hop for all particles in saltation, the average particle streamwise acceleration, 〈*a*_{p}〉, is approximated by $\langle \mathrm{\Delta}{v}_{x,\mathrm{s}}\rangle /\langle {t}_{\mathrm{p}}\rangle $. In this way, the transport rate in Eq. (1) can be expressed by

According to the balance of linear momentum applied to the particles in saltation described in Sect. 2.2, the ratio between the total aerodynamic force along the streamwise direction applied to all particles in saltation and their weight is given by $\langle {a}_{\mathrm{p}}\rangle /g$. Bagnold (1973) argued that this ratio between the applied horizontal force and the weight of the mass in motion represents a mean friction coefficient. This friction coefficient is expected to be of a similar nature to the coefficient of static friction obtained between two solid surfaces in contact. Taking into account that the real (microscopic) contact between two surfaces is actually discontinuous, Bagnold (1973) claims that this analogy holds in the case of saltating particles, in which the contact between the particles and the surface is intermittent. Therefore, similar to a coefficient of static friction, the ratio $\langle {a}_{\mathrm{p}}\rangle /g$ is expected to be a function of the particle–bed interaction only and, consequently, invariant with respect to *u*_{∗}.

In the saltation model proposed by Pomeroy and Gray (1990), an expression for the transport rate of saltating snow is also obtained by estimating the ratio $\langle {a}_{\mathrm{p}}\rangle /g$. However, in contrast to Bagnold (1973), 〈*a*_{p}〉 was found to increase linearly with *u*_{∗}. In addition, 〈*a*_{p}〉 was found to be highly variable with the snow surface characteristics. In the model, 〈*a*_{p}〉 was not estimated from the direct measurements of 〈Δ*v*_{x,s}〉 and 〈*t*_{p}〉, but it was deduced from fitting experimental measurements of the particle mass flux in saltation with Eq. (2). Therefore, the resultant scaling of 〈*a*_{p}〉 with *u*_{∗} is a function of the parameterizations considered for *τ*_{s} and 〈*v*_{p}〉, as well as of the accuracy of the mass flux measurements and the respective vertically integrated value, *Q*.

The mass flux measurements considered in the model were performed at Saskatoon, Canada, using an optoelectronic snow particle detector placed at a height of approximately 2 cm above the snow surface. This device counts the number of particles crossing an infrared laser beam. The mass flux is then estimated by assuming a particle size distribution given by a two-parameter gamma distribution, as suggested by Schmidt (1981), and a particle density equal to the density of ice (∼ 917 kg m^{−3}). The mass flux in saltation is assumed to be uniform in height. Therefore, the transport rate is estimated from the product of the measured mass flux and the height of the saltation layer, *h*_{salt} [m]. As addressed in Sect. 2.3.3, the vertical profile of the particle mass flux exhibits strong gradients close to the surface. In this way, the saltation mass flux is hardly defined by a single-point measurement and the assumption of a constant mass flux along the saltation layer might lead to errors in the assessment of the transport rate.

In summary, Pomeroy and Gray (1990) arrived at the conclusion that 〈*a*_{p}〉 scales with *u*_{∗} by assuming that the mass flux is uniform in height, that *h*_{salt} scales with ${u}_{\ast}^{\mathrm{2}}$ (Owen, 1964) and that 〈*v*_{p}〉 is invariant with respect to *u*_{∗}. As a consequence, they concluded that *Q* scales linearly with *u*_{∗}. This result is not in agreement with most saltation models and was shown to highly underestimate the transport rate in comparison to other expressions proposed in the literature (Melo et al., 2022). Notwithstanding, this model is still widely used, specifically in the drifting snow schemes of RACMO (van Wessem et al., 2018) and MAR (Amory et al., 2021). The model of Pomeroy and Gray (1990) is one of the few parameterizations for the transport rate that take into account field measurements of snow saltation. However, it is highly dependent on the assumptions made for the particle mass flux and streamwise velocity profiles. These will be discussed in Sect. 2.3.

If both 〈Δ*v*_{x,s}〉 and 〈*t*_{p}〉 are invariant with respect to *u*_{∗}, as suggested by several models and measurements (see Sect. 2.2.2 and 2.2.3), 〈*a*_{p}〉 is mainly expected to be invariant with respect to *u*_{∗}, as hypothesized earlier by Bagnold (1973). Taking into account that the direct measurement of 〈*a*_{p}〉 is challenging (e.g., Araoka and Maeno, 1981), the variation in *Q* with *u*_{∗} and the snow surface characteristics might be better understood from the direct study of 〈Δ*v*_{x,s}〉, 〈*t*_{p}〉 and 〈*v*_{p}〉.

## 2.3 The vertical profiles of particle concentration, streamwise velocity and mass flux

The particle mass flux in saltation, *q* [$\mathrm{kg}\phantom{\rule{0.25em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.25em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$], corresponds to the mass flow rate of particles moving in saltation per unit cross section area. The particle concentration, *c* [kg m^{−3}], is the mass of particles in saltation per unit volume. By definition, these two quantities are related by $q=c{\stackrel{\mathrm{\u203e}}{v}}_{x}$, where ${\stackrel{\mathrm{\u203e}}{v}}_{x}$ [m s^{−1}] is the volume-averaged particle streamwise velocity (henceforth, particle streamwise velocity) at which the respective mass of particles in saltation travels across a given cross section area (note that the average of ${\stackrel{\mathrm{\u203e}}{v}}_{x}$ along the saltation layer equals the average particle streamwise velocity, 〈*v*_{p}〉, defined in Sect. 2.2). During steady-state saltation, *q*, *c* and ${\stackrel{\mathrm{\u203e}}{v}}_{x}$ are invariant along the horizontal directions, but they are expected to vary with the height above the surface. In the following subsections, we discuss the parameterizations used to model these quantities and the experimental, theoretical and numerical results that contributed to their development.

### 2.3.1 The particle concentration

An analytic expression for the vertical profile of particle concentration was derived by Kawamura (1951). In the model, the saltation dynamics is simplified by neglecting the aerodynamic forces on the vertical motion of saltating particles, by assuming that the probability distribution of the vertical component of the particle ejection velocity follows a half-normal distribution and that the saltating particles are uniformly sized. In this way, from a mass balance between the surface and a given height *z*, the particle concentration at each height can be estimated as a function of the probability distribution of the vertical ejection velocity and the vertical mass flux of particles leaving the surface. In Gordon et al. (2009), it is shown that the analytic expression of Kawamura (1951) is well approximated by an exponential decay of the form

where *c*_{o} [kg m^{−3}] is the particle concentration at the surface (by extrapolating the exponential profile down to the surface) and *δ*_{c} [m] is the decay height of the particle concentration profile. The decay height, *δ*_{c}, defines the height above the surface for which the particle concentration is reduced to ${c}_{\mathrm{o}}/\mathrm{exp}\left(\mathrm{1}\right)$. The theoretical expression proposed by Kawamura (1951) tends to infinity at the surface. In contrast, Eq. (3) is expected to underestimate the particle concentration near the surface and to provide a good approximation only for *z*>〈*h*_{p}〉, where 〈*h*_{p}〉 [m] is the average hop height of particles in saltation.

The direct measurement of the particle concentration is challenging because it requires the measurement of the number and size of saltating particles in a well-defined volume (Creyssels et al., 2009). Gordon et al. (2009) measured the number density of snow particles in saltation in the first 9 cm above the surface over sea ice in Franklin Bay, Canada. Images of saltating snow were recorded with a camera, which captured the light reflected by snow particles traveling in front of a black background. The vertical profile of particle concentration was estimated by specifying the width over which the particles were detected (volume depth), by assuming that the particle size distribution is given by a two-parameter gamma distribution and that the saltating particles have the density of ice. The resultant particle concentration profiles, obtained for friction velocities ranging from 0.25 to 0.5 m s^{−1} and different snow surface conditions, followed mainly an exponential decay as predicted by the theoretical model of Kawamura (1951). However, when the particle concentration was low, the data were more scattered, and deviations from the exponential decay were seen above a height of approximately 2 cm above the surface.

From fitting the measured vertical profiles with Eq. (3), the respective decay height was estimated. According to Gordon et al. (2009), *δ*_{c} is expected to scale with the average hop height, 〈*h*_{p}〉. In the absence of aerodynamic forces, the hop height of a particle in saltation scales with ${{v}_{z}^{\mathrm{e}}}^{\mathrm{2}}/g$. Therefore, *δ*_{c} is expected to scale with $\langle {v}_{z}^{\mathrm{e}}{\rangle}^{\mathrm{2}}/g$. Owen (1964) suggested that $\langle {v}_{z}^{\mathrm{e}}\rangle $ scales with *u*_{∗} and, therefore, that 〈*h*_{p}〉 scales with ${u}_{\ast}^{\mathrm{2}}$. However, the measurements of Gordon et al. (2009) revealed a weak correlation between *δ*_{c} and ${u}_{\ast}^{\mathrm{2}}$. Indeed, for each set of measurements conducted on the same day (similar snow surface conditions), the decay height did not increase monotonically with the friction velocity.

Similar results were found by Creyssels et al. (2009) after performing wind tunnel experiments with uniformly sized sand particles. As mentioned in Sect. 2.2.2, the experiments were performed with a high-speed camera illuminated by a laser sheet. The particle concentration was computed in the first 5 cm above the surface for friction velocities ranging from 0.24 to 0.67 m s^{−1}. In the calculations, the volume depth was assumed equal to the width of the laser. The vertical profiles of particle concentration were found to follow an exponential decay of the form of Eq. (3), and the respective decay height was found to be invariant with respect to the friction velocity. The conclusions of Gordon et al. (2009) and Creyssels et al. (2009) regarding the scaling of *δ*_{c} with *u*_{∗} are in agreement with the assumption that the ejection velocity is invariant with respect to *u*_{∗}, which is in line with most numerical models and experimental measurements developed over erodible surfaces (see Sect. 2.2.2).

Ho et al. (2011) used a similar technique to measure the particle concentration in the first 6 cm above the surface in a wind tunnel. In their experiments, both a bed composed of loose sand grains and a rigid (non-erodible) surface were used. In agreement with the experiments of Creyssels et al. (2009), the particle concentration over the erodible surface was found to decrease with increasing height following an exponential decay, and the computed decay height was invariant with respect to *u*_{∗}. In contrast, when saltation was allowed to develop over the rigid surface, the decay height was found to be higher and to increase linearly with *u*_{∗}. This result is related to the higher values of the particle slip velocity obtained over the rigid surface in comparison to the erodible one and to the increase in the particle slip velocity with the increase in *u*_{∗} obtained over the rigid surface, mentioned in Sect. 2.2.2.

Taking into account the challenges and uncertainties in the particle concentration measurements, especially in the field, measurements of particle mass flux are more frequent (e.g., Namikas, 2003; Martin and Kok, 2017; Nishimura and Nemoto, 2005). Therefore, the general approach in drifting snow schemes is to derive the particle concentration as a function of the particle streamwise velocity and mass flux (Amory et al., 2021; van Wessem et al., 2018; Vionnet et al., 2014; Sharma et al., 2023). These two quantities are discussed in the next sections.

### 2.3.2 The particle streamwise velocity

The particle streamwise velocity of saltating sand was measured by Creyssels et al. (2009) and Ho et al. (2011) in the first 5 and 3 cm above the surface, respectively (see details in Sect. 2.2.2). In this region of the saltation layer, the particle streamwise velocity, ${\stackrel{\mathrm{\u203e}}{v}}_{x}$, was found to increase linearly with height. Over a bed of loose sand grains, the effect of the friction velocity on the particle slip velocity was almost negligible, as discussed in Sect. 2.2.2. Above the surface, ${\stackrel{\mathrm{\u203e}}{v}}_{x}$ was found to increase slightly with the rise in *u*_{∗}. Aksamit and Pomeroy (2016), also using the PTV technique, measured the velocity of saltating snow in the first 2.5 cm above the surface in a snow-covered area. A linear increase in ${\stackrel{\mathrm{\u203e}}{v}}_{x}$ with height was seen in the first 1 cm above the surface. In addition, ${\stackrel{\mathrm{\u203e}}{v}}_{x}$ was found to increase slightly with the friction velocity.

The particle streamwise velocity can also be measured with a snow particle counter (SPC) (Nishimura et al., 2014). An SPC is an optical device that detects the particles that cross a 25 mm long, 2 mm high and 0.5 mm wide laser beam. In general, this device is used to retrieve the number and average diameter of all particles crossing the laser beam during a defined time period. Therefore, it is mainly used to estimate the particle mass flux. The particle streamwise velocity can be computed by evaluating the time that each particle takes to cross the laser beam. Nishimura et al. (2014) measured the streamwise velocity of saltating snow along the first 10 cm above the surface in a wind tunnel for friction velocities ranging from 0.37 to 0.63 m s^{−1} using this technique. In agreement with the previous studies, the particle streamwise velocity was found to increase as the height above the surface and the friction velocity increase. However, in contrast to the previous PTV measurements, most of the measurement points were located above a height of 5 cm. In this way, the increase in ${\stackrel{\mathrm{\u203e}}{v}}_{x}$ with increasing height did not follow a linear trend but rather a logarithmic one. Similar particle streamwise velocity profiles were found in the wind tunnel measurements performed by Nishimura and Hunt (2000), who recorded saltating snow particles and ice spheres in the first 10 cm above the surface with a video system, and by Rasmussen and Sørensen (2008), who used laser Doppler velocimetry (LDV) to measure the streamwise velocity of saltating sand in the first 8 cm above the surface.

Numerical models that represent the surface processes of splash and rebound have also obtained an increase in the particle streamwise velocity with height (e.g., Kok and Renno, 2009; Durán et al., 2011). As expected, in the first 1 cm above the surface, the particle velocity is mainly invariant with respect to *u*_{∗}, and, at higher elevations, the particle streamwise velocity increases with the increase in *u*_{∗}.

In the model of Pomeroy and Gray (1990), the particle streamwise velocity is assumed constant in height and is characterized by its average value 〈*v*_{p}〉. 〈*v*_{p}〉 is assumed invariant with respect to *u*_{∗} and proportional to the impact threshold friction velocity, ${u}_{\ast ,\text{it}}$. This assumption can be found in the drifting snow schemes of RACMO, MAR, Meso-NH and CRYOWRF to compute the particle concentration at the top of the saltation layer as a function of the particle streamwise velocity and mass flux (Lenaerts et al., 2012; Amory et al., 2021; Vionnet et al., 2014; Sharma et al., 2023). The assumption that 〈*v*_{p}〉 scales with ${u}_{\ast ,\text{it}}$ is justified by the observations of Bagnold (1941) regarding the existence of a focus point in the wind velocity profile. According to his observations, during steady-state saltation, the wind velocity inside the saltation layer obtained for different values of *u*_{∗} converge to the same value at a given height close to the surface: the focus point. The wind velocity at the focus point scales with the impact threshold friction velocity and is therefore solely dependent on the bed characteristics. In this way, the particle streamwise velocity in the vicinity of the focus point is also expected to scale with ${u}_{\ast ,\text{it}}$, which according to Owen (1964) is the value of the surface friction velocity during saltation (see Sect. 2.2.1). The existence of a focus point and the fact that the particle streamwise velocity close to the surface is invariant with respect to *u*_{∗} are known features of the saltation system, obtained by different models that represent the surface processes of rebound and splash (e.g., Kok and Renno, 2009; Durán et al., 2011; Melo et al., 2022). However, as previously discussed, experimental measurements performed with both sand and snow show that the particle streamwise velocity above a height of approximately 3 cm increases with the increase in *u*_{∗}. In this way, the parameterization used for ${\stackrel{\mathrm{\u203e}}{v}}_{x}$, in particular, its scaling with ${u}_{\ast ,\text{it}}$ and *u*_{∗}, must be defined according to the height at which the particle streamwise velocity (and, subsequently, the particle concentration) is being computed.

The snow surface characteristics are also expected to influence the particle streamwise velocity, as previously discussed in Sect. 2.2.2. For instance, for the same range of friction velocities, Aksamit and Pomeroy (2016) found a higher particle streamwise velocity in the first 2.5 cm above the surface when saltation developed over a surface of fresh snow than over surfaces of old snow or fine decomposing grains. One can hypothesize that this is related to the smaller density and dendritic shape of the particles in saltation. The particles ejected from a fresh snow surface are not fully rounded and will continuously fragment and get denser as the particles collide with the bed (Comola et al., 2017). Other snow surface characteristics, like the particle diameter and the existence of cohesive bonds between the snow particles, can also influence the saltation dynamics. When modeling the impact and splash of saltating snow particles, Comola and Lehning (2017) took into account the effect of inter-particle cohesion. In their model, a fraction of the kinetic energy of the impacting particle is used to break the cohesive bonds between each ejected particle and the surrounding ones. Melo et al. (2022) used this splash model to simulate saltating snow and assessed the effect of the particle diameter and inter-particle cohesion on the particle streamwise velocity. Considering a bed of uniformly sized particles, the particle streamwise velocity was found to decrease with the increase in the particle diameter. Moreover, by increasing the energy required to break the cohesive bonds, a decrease in the particle concentration and an increase in the particle streamwise velocity were found in the first 5 cm above the surface. In the model, the existence of inter-particle bonds leads to a decrease in the number of particles aloft. This small number of particles is able to extract more momentum from the flow and to reach higher velocities. The same trend was found by Comola et al. (2019a) when studying the effect of inter-particle cohesion on steady-state saltation, using the discrete element method (DEM) and a simple description of the wind velocity field.

### 2.3.3 The particle mass flux

The vertical profile of particle mass flux is expected to follow an exponential decay of the same form of Eq. (3) (e.g., Nalpanis et al., 1993; Namikas, 2003; Martin and Kok, 2017):

where *q*_{o} [$\mathrm{kg}\phantom{\rule{0.25em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.25em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$] is the surface mass flux (by extrapolating the exponential profile down to the surface) and *δ*_{q} [m] is the decay height of the particle mass flux profile. By definition, the integral of *q* along the height yields the transport rate, *Q*. Therefore, it follows that ${q}_{\mathrm{o}}=Q/{\mathit{\delta}}_{\mathrm{q}}$.

Experiments developed with sand show that the mass flux profile follows an exponential decay at moderate friction velocities but does not at low values of *u*_{∗}, close to the fluid threshold friction velocity (Nalpanis et al., 1993; Namikas, 2003). In addition, it was shown that the exponential decay tends to underestimate the mass flux in the first 2 cm above the surface (Namikas, 2003; Bauer and Davidson-Arnott, 2014). This might pose a problem when relating the fitted values of *q*_{o} and *δ*_{q} to the transport rate, *Q*. Experiments developed with snow, both in wind tunnels (Sugiura et al., 1998; Sato et al., 2001) and in the field (Takeuchi, 1980; Nishimura and Nemoto, 2005), revealed an exponential decay of the mass flux profile in approximately the first 10 cm above the surface. Near the top of the saltation layer, deviations are seen due to the saltation–suspension transition (Gordon et al., 2009). The exponential decay of the particle mass flux with height contrasts with the uniform profile assumed by Pomeroy and Gray (1990), as discussed in Sect. 2.2.4.

Similarly to the decay height presented in Eq. (3), *δ*_{q} is also thought to scale with the average particle hop height in saltation and, therefore, with $\langle {v}_{z}^{\mathrm{e}}{\rangle}^{\mathrm{2}}/g$ (Nishimura and Hunt, 2000). Based on the assumption made by Owen (1964) regarding the scaling of $\langle {v}_{z}^{\mathrm{e}}\rangle $ with *u*_{∗}, Nishimura and Hunt (2000) assumed that *δ*_{q} must scale with ${u}_{\ast}^{\mathrm{2}}/g$. This parameterization for *δ*_{q} is currently used in Meso-NH (Vionnet et al., 2014) and CRYOWRF (Sharma et al., 2023) to compute the particle mass flux in saltation.

The decay height is estimated from fitting Eq. (4) to measured mass flux profiles. In Fig. 2, the values of *δ*_{q} obtained from different experiments developed with sand and snow are presented as a function of the friction velocity, *u*_{∗} (Sugiura et al., 1998; Namikas, 2003; Nishimura and Nemoto, 2005; Martin and Kok, 2017). The expression proposed by Nishimura and Hunt (2000) is also presented for comparison. The measurements of Sugiura et al. (1998) were performed in a wind tunnel using preserved natural snow, which was disintegrated into individual particles prior to the experiments. The particle mass flux was measured between 1.6 and 6.1 cm above the surface using an SPC. By detecting the number and average diameter of saltating particles, the mass flux is computed by approximating the snow particles to spheres and by assuming a given particle density (generally, 917 kg m^{−3}). As seen in Fig. 2, the values of *δ*_{q} obtained for friction velocities ranging from 0.15 to 0.39 m s^{−1} were found to increase with ${u}_{\ast}^{\mathrm{2}}$, as suggested by Nishimura and Hunt (2000).

The results of Sugiura et al. (1998) are in agreement with the wind tunnel experiments of Sato et al. (2001) performed with two different snow covers. Similarly to the experiments of Sugiura et al. (1998), snow particles were obtained by disintegrating a layer of preserved natural snow. Experiments were performed immediately after spreading the snow particles over the wind tunnel floor (loose snow cover) and after some time, which allowed sintering to take place (hard snow cover). The vertical profile of particle mass flux followed an exponential decay over both snow surfaces, and the decay height was found to increase with the wind velocity. However, significantly different values of *δ*_{q} were obtained with each snow cover: for the same range of wind velocities, *δ*_{q} varied between approximately 0.9 and 3 cm over the loose snow cover and between approximately 3 and 8 cm over the hard snow cover. Even though the hard snow cover cannot be considered an erodible surface (particle ejection from the surface was rare, and the mass flux in saltation was highly dependent on the seeding rate of snow particles at the upwind end of the tunnel), the results of Sato et al. (2001) highlight the influence of the snow surface characteristics on the saltation dynamics and, in particular, on the decay height. The decay height obtained by Sato et al. (2001) is not presented in Fig. 2 because the respective value of *u*_{∗} was not provided by the authors.

The mass flux measurements of Nishimura and Nemoto (2005) were performed at Mizuho Station, Antarctica, using SPCs. They reported mass flux profiles obtained for friction velocities ranging from 0.21 to 0.56 m s^{−1} (please note that the mass flux measurements reported by Nishimura and Nemoto (2005) in Fig. 6 must be multiplied by 10 to be in agreement with the units; Kouichi Nishimura, personal communication, 22 November 2023). The respective *δ*_{q} values are not given in the article, but these can be estimated by fitting the mass flux measurements between heights of 2 and 12 cm to Eq. (4). From Fig. 2, one can conclude that *δ*_{q} increases with *u*_{∗}. In addition, considering the three lowest friction velocities, a fair agreement is obtained between the computed values of *δ*_{q} and the expression of Nishimura and Hunt (2000). At higher friction velocities, higher values of *δ*_{q} are obtained in comparison with the quadratic expression. However, without a detailed characterization of the snow cover during the different events (which spanned a time period of 1.5 months), it is not possible to evaluate the effect of the friction velocity alone.

When performing field and wind tunnel experiments with sand, a different trend is obtained. For instance, Namikas (2003) measured the particle mass flux along the first 35 cm above the surface using wedge-shaped sediment collectors. The measurements were performed for friction velocities ranging from 0.27 to 0.63 m s^{−1}. From fitting the measured profiles with Eq. (4), similar values of *δ*_{q} were obtained. Therefore, it was concluded that *δ*_{q} is invariant with respect to *u*_{∗}. The same conclusion was reached by Martin and Kok (2017) after performing field experiments in three different locations: Oceano, Rancho Guadalupe and Jericoacoara. The particle mass flux was measured with Big Spring Number Eight (BSNE) sand traps between heights of approximately 15 and 60 cm above the surface. Even though this height range is above the saltation layer in the case of snow saltation, the mass flux of saltating sand still follows an exponential decay in this region (Martin et al., 2018). Based on their results, Martin and Kok (2017) suggest that *δ*_{q} rather scales with the particle diameter. This would justify the agreement between the measurements of Martin and Kok (2017) at Oceano and those of Namikas (2003), which were performed at a similar location. Bigger sand grains were reported at Rancho Guadalupe and Jericoacoara in comparison to Oceano (Martin and Kok, 2017). Similarly, Nalpanis et al. (1993) could not verify the scaling of *δ*_{q} with ${u}_{\ast}^{\mathrm{2}}$ after performing wind tunnel experiments with sand and numerical simulations.

By measuring the particle concentration and streamwise velocity, Creyssels et al. (2009) found that the mass flux evolution is mainly driven by the particle concentration profile. Therefore, even though the particle streamwise velocity was found to increase linearly with height in the first 5 cm above the surface, *δ*_{c} and *δ*_{q} were found to be approximately equal. In this way, the scaling of *δ*_{c} with *u*_{∗}, reported in Sect. 2.3.1, is expected to reflect the scaling of *δ*_{q} with *u*_{∗}. In agreement with the mass flux measurements performed with sand, Creyssels et al. (2009) and Gordon et al. (2009) did not arrive at an increase in *δ*_{c} with the increase in *u*_{∗}. In contrast, when saltation was allowed to develop over a rigid surface, Ho et al. (2011) obtained a linear increase in *δ*_{c} with *u*_{∗}. Despite the findings of Creyssels et al. (2009) regarding the equality between *δ*_{c} and *δ*_{q}, the values of *δ*_{q} presented in Fig. 2 and the values of *δ*_{c} presented in the works of Creyssels et al. (2009), Ho et al. (2011) and Gordon et al. (2009) differ by approximately 1 order of magnitude. In this way, the comparison between these two quantities must be performed with care.

The scaling of *δ*_{q} with *u*_{∗} can also be evaluated from vertical profiles of particle number flux. Considering saltating particles of uniform size, the number flux is proportional to the mass flux and, therefore, is expected to follow an exponential decay of the form of Eq. (4). Aksamit and Pomeroy (2016) computed the number flux of snow particles in saltation in the first 1 cm above the surface from the recordings of high-speed images of saltating snow in a snow-covered area. Even though the observed particle size distribution was not uniform but characterized by a gamma distribution, the decay height of the number flux profile is not expected to deviate significantly from *δ*_{q}. Similarly to the results obtained by Gordon et al. (2009) over sea ice, a clear scaling of the decay height with *u*_{∗} was not found: it either increased or decreased with the rise in the friction velocity depending on the snow cover characteristics (please note that the values reported by Aksamit and Pomeroy (2016) in Fig. 7d–f do not correspond to *l*_{v} [mm], named decay length, but to $\mathrm{1000}/{l}_{\mathrm{v}}$; Nikolas O. Aksamit, personal communication, 11 February 2023). In addition, the values obtained are 1 order of magnitude lower than those reported by Sugiura et al. (1998) for the same range of friction velocities and of the same order of magnitude as those reported in Gordon et al. (2009).

The uncertainty regarding the scaling of *δ*_{q} with *u*_{∗} can, in part, be linked to the uncertainty regarding the scaling of $\langle {v}_{z}^{\mathrm{e}}\rangle $ with *u*_{∗} when saltation develops over a snow surface (see Sect. 2.2.2). The idea that *δ*_{q} is invariant with respect to *u*_{∗} is in agreement with the assumption that the ejection velocity is invariant with respect to *u*_{∗}. The increase in *δ*_{q} with the rise in *u*_{∗} obtained by Sugiura et al. (1998) and Nishimura and Nemoto (2005) is not in agreement with this assumption and contrasts with the results of Gordon et al. (2009) and Aksamit and Pomeroy (2016) regarding the scaling of the decay height of the particle concentration and number flux profiles with respect to *u*_{∗}. These contradictory results remain to be understood. In addition, more experimental studies are needed to understand and quantify the effect of the snow surface characteristics on *δ*_{q}.

In this section, we use a numerical model to investigate the main quantities of interest in snow saltation modeling, previously discussed in Sect. 2. The numerical model comprises a large-eddy simulation flow solver; a Lagrangian stochastic model for the particle trajectories; and a set of parameterizations that represent the processes of aerodynamic entrainment, rebound and splash (Sharma et al., 2018; Comola et al., 2019b; Sigmund et al., 2022; Melo et al., 2022). A brief description of the model is given in Sect. 3.1. Detailed numerical models, like the one described, cannot be used in atmospheric models, which require simple and computationally light algorithms. Nevertheless, they provide an excellent test bench for simple parameterizations and their underlying assumptions. The description of the numerical setup is presented in Sect. 3.2, and the main results are presented and discussed in Sect. 3.3.

## 3.1 Model description

The turbulent wind field is described by the continuity and Navier–Stokes equations and is assumed incompressible. Following the large-eddy simulation (LES) technique, a filtering operator is applied to the equations in order to resolve the turbulent structures larger than the computational grid size. The turbulent structures smaller than the grid size are not resolved, but their effect on the larger structures is parameterized with a subgrid-scale (SGS) model. The numerical implementation of the LES model is based on the work of Albertson and Parlange (1999), which was continuously developed for different boundary layer studies (e.g., Giometto et al., 2016, 2017; Sharma et al., 2017). In the model, we consider the Lagrangian-averaged scale-dependent SGS model proposed by Bou-Zeid et al. (2005). The computational domain is of cuboid shape, in which periodic boundary conditions are imposed at the vertical walls. Impermeability and zero vertical gradients are assumed at the top boundary. At the bottom boundary, the logarithmic law of the wall and the impermeability condition are specified. The flow is driven by a constant streamwise pressure gradient equal to $-{\mathit{\rho}}_{\mathrm{a}}{u}_{\ast}^{\mathrm{2}}/{L}_{z}$, where *L*_{z} [m] is the domain height. In the absence of saltating particles, this imposed pressure gradient guarantees a time-averaged surface shear stress equal to ${\mathit{\rho}}_{\mathrm{a}}{u}_{\ast}^{\mathrm{2}}$.

The particle trajectories within the turbulent flow are modeled with a Lagrangian stochastic model (LSM) (Comola, 2017). The particle motion is computed with Newton's second law, considering the effects of gravity and aerodynamic drag. The particles are assumed spherical, and the drag coefficient is computed with the expression of Schiller and Nauman (1933). Other forces applied to the particles (e.g., aerodynamic lift, electrostatic force and virtual mass) are neglected in the model. The feedback of the particles on the flow momentum is taken into account by the addition of a reaction force in the filtered Navier–Stokes equations. It is equal in magnitude and opposite to the aerodynamic drag applied to the particles.

The surface processes of aerodynamic entrainment, rebound and splash occur at the bottom boundary, which is assumed erodible and composed of particles with a predefined size distribution. The surface processes are described by physical and statistical models, as proposed by Groot Zwaaftink et al. (2014) and Comola and Lehning (2017). For instance, aerodynamic entrainment of surface particles is expected to occur when the surface shear stress is greater than the fluid threshold (Bagnold, 1941):

where *A*=0.1 is the fluid threshold coefficient, *ρ*_{p} [kg m^{−3}] is the particle density and $\stackrel{\mathrm{\u203e}}{d}$ [m] is the average particle diameter of the erodible bed. The number of entrained particles is proportional to the difference *τ*_{s}−*τ*_{ft}
(Anderson and Haff, 1991; Doorschot and Lehning, 2002). Following Clifton and Lehning (2008), it is assumed that the ejection velocity of an aerodynamically entrained particle varies according to a lognormal distribution. The mean and standard deviation of the distribution are assumed to be proportional to ${u}_{\ast ,\mathrm{s}}$. The initial particle concentration at the surface is considered high enough so that there is never a shortage in the supply of erodible particles.

Following a particle impact with the granular surface, the impacting particle can rebound, and some particles can be splashed from the surface. The probability of rebound increases as the particle impact velocity increases (Anderson and Haff, 1991), and the rebound velocity is given by $\sqrt{{\mathit{\u03f5}}_{\mathrm{r}}}\left|{\mathit{v}}^{\mathrm{i}}\right|$, where *ϵ*_{r}=0.25 is the fraction of kinetic energy retained by the rebounding particle (restitution coefficient) and |*v*^{i}| [m s^{−1}] is the modulus of the impact velocity (Doorschot and Lehning, 2002). The number of splashed particles, *N*, is computed with the model of Comola and Lehning (2017). In their model, the balances of kinetic energy and momentum of the collision are resolved by statistically representing the kinetic energy and momentum of the splashed particles by their mean values. As a result, the number of splashed particles varies with the characteristics of the impacting particle (mass and velocity) and the statistical characteristics of the splashed particles (mean and standard deviation of the particle size and ejection velocity). We assume that the ejection velocity of splashed particles varies according to an exponential distribution. Its mean value is a function of the particle impact velocity and is given by $\mathrm{0.25}{\left|{\mathit{v}}^{\mathrm{i}}\right|}^{\mathrm{0.3}}$ (Sharma et al., 2018). In the model, *N* is also a function of other parameters that describe the granular bed as, for example, the restitution coefficient; the energy required to break the cohesive bonds between each ejected particle and the surrounding particles; and the fraction of impact kinetic energy and momentum that is lost due to friction. The full description of the LES-LSM model can be found in previous works (e.g., Sharma et al., 2018; Comola et al., 2019b; Melo et al., 2022).

## 3.2 Numerical settings

The computational domain consists of a cube with a 6.4 m side length, uniformly discretized in each horizontal direction in 64 cells. The vertical direction is discretized in 128 cells following a hyperbolic stretching, which guarantees a more refined mesh close to the bottom boundary. The simulations are performed for 350 s. The flow is allowed to develop over the first 25 s of the simulations. The start of surface erosion is only allowed after this first time period. The time step of the LES solver is set to 5 × 10^{−5} s.

In order to reproduce the snow surface characteristics reported in Sugiura et al. (1998), the particle size distribution of the granular bed is assumed normal, with a mean diameter of 360 µm and a standard deviation of 140 µm. A minimum and maximum particle diameter is imposed, equal to 30 µm and 2 mm, respectively. In the experiments of Sugiura et al. (1998), the wind tunnel was maintained at −15 °C. Therefore, in the simulations, the air density and kinematic viscosity, *ν* [m^{2} s^{−1}], are set to the respective values at that temperature (*ρ*_{a} = 1.37 kg m^{−3}, *ν* = 1.2 × 10^{−5} m^{2} s^{−1}). Moreover, snow sublimation is neglected. Taking into account that the snow layer was composed of preserved natural snow that was disintegrated into individual particles prior to the experiments, the effect of inter-particle cohesion is neglected in the simulations. The remaining rebound and splash parameters are set to the values reported in Melo et al. (2022). In addition, the particle density is assumed equal to 917 kg m^{−3}, as suggested by Sugiura et al. (1998). The roughness length, *z*_{o}, is not reported in Sugiura et al. (1998), and a value of 10^{−4} m is assumed (Clifton et al., 2006).

A set of seven simulations is performed for which the friction velocity is varied between 0.15 and 0.7 m s^{−1}. In order to decrease the computational time, the time step of the LSM solver is set to 10^{−4} s (twice the LES time step) for the friction velocities ranging between 0.15 and 0.5 m s^{−1} and to 2 × 10^{−4} s (4 times the LES time step) for the higher friction velocities. In addition, the particles in saltation are grouped in parcels, composed of particles with equal diameter and following the same trajectory. Particles from the same parcel were aerodynamically entrained or splashed at the same surface location and time step. In order to guarantee that the simulation results are not affected by this approximation, the number of particles per parcel is set to 5, 50 and 100 for the friction velocities equal to 0.15, 0.23 and 0.3 m s^{−1}, respectively, and to 200 for the higher friction velocities.

## 3.3 Results and discussion

The vertical profiles of particle mass flux and velocity are computed by dividing the computational domain in horizontal layers. At each layer, the mass flux is computed by

where *z*_{k} [m] and Δ*z*_{k} [m] are the average height above the surface and the thickness of layer *k*, *N*_{k} is the total number of particles in layer *k*, *m*_{n} [kg] and ${v}_{{x}_{n}}$ [m s^{−1}] are the mass and the instantaneous streamwise velocity of the *n*th particle in the respective layer, and *L*_{x} [m] and *L*_{y} [m] are the domain length and width. The particle streamwise velocity at each layer is given by the average instantaneous particle velocity considering all particles in the layer. The vertical profiles are time-averaged over the last 100 s of the simulations. During this time period, negligible changes in the total mass of particles aloft were obtained, and saltation was assumed to be in steady state.

The obtained vertical profiles of particle mass flux are presented in Fig. 3, together with the SPC measurements performed by Sugiura et al. (1998). The results obtained for friction velocities ranging from 0.15 to 0.5 m s^{−1} are presented in a linear plot (Fig. 3a), and those obtained for all simulated friction velocities are presented in a semi-logarithmic plot (Fig. 3b). Curves fitting the mass flux profiles with Eq. (4) are also presented. It can be seen that most of the mass flux profiles follow an exponential decay in approximately the first 8 cm above the surface. The only exception is the mass flux profile obtained with the lowest friction velocity of 0.15 m s^{−1}. In this case, the friction velocity outside the saltation layer is slightly lower than the assumed fluid threshold friction velocity (from Eq. (5), ${u}_{\ast ,\text{ft}}$ = 0.154 m s^{−1}). Therefore, saltation is highly intermittent. Sugiura et al. (1998) have also reported a fluid threshold friction velocity close to 0.15 m s^{−1}. During the experiments, saltation was only sustained at this low friction velocity with the addition of seeding particles at the upwind end of the wind tunnel. For the remaining friction velocities, deviations are seen between the modeled mass flux profiles and the exponential decay in approximately the first 1 cm above the surface and above heights of 6–12 cm, depending on the friction velocity. The deviations from an exponential decay above heights of 6–12 cm are due to the prevalence of particles in suspension. As discussed in Sect. 2.3.3, these trends were also observed in field measurements of snow and sand saltation (Namikas, 2003; Bauer and Davidson-Arnott, 2014; Nishimura and Nemoto, 2005).

When comparing the modeled mass flux profiles with the SPC measurements of Sugiura et al. (1998) obtained at the same friction velocities, it is clear that the numerical model is significantly underestimating the mass flux in saltation. This discrepancy could not be solved by decreasing the roughness length from 10^{−4} to 10^{−5} m, which is more in agreement with previous wind tunnel experiments developed with snow (Nishimura et al., 2014). Negligible improvements were also obtained by decreasing the domain height to half the height of the wind tunnel (please note that the top boundary of the domain is a symmetry plane). Despite the qualitative agreement, a quantitative agreement between simulation results and wind tunnel measurements certainly requires an adjustment in the parameters that control the surface processes of aerodynamic entrainment, rebound and splash. Even though a fitting exercise could be done to achieve a better agreement, we believe that detailed measurements of the near-surface particle–bed interaction and of the snow particle characteristics in saltation are needed for a further understanding of the near-surface processes and the development of improved parameterizations.

From fitting the computed mass flux profiles to Eq. (4), the decay height can be obtained. It is presented in Fig. 4 as a function of the friction velocity. The values obtained by Sugiura et al. (1998) after fitting the measured mass flux profiles to Eq. (4) and the expression proposed by Nishimura and Hunt (2000) are also presented for comparison. For friction velocities ranging from 0.15 to 0.3 m s^{−1}, there is a fair agreement between the decay height obtained from measured and modeled mass flux profiles. The decay height increases with the increase in the friction velocity and does not deviate considerably from the evolution proposed by Nishimura and Hunt (2000). However, for friction velocities greater than 0.3 m s^{−1}, the values of decay height obtained from the modeled mass flux profiles increase with the increase in the friction velocity at a much lower rate. In particular, the value obtained for ${u}_{\ast}=\mathrm{0.39}\phantom{\rule{0.25em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ is lower than the value obtained by Sugiura et al. (1998), and a significant deviation between the modeled values and the expression proposed by Nishimura and Hunt (2000) is seen.

Overall, the decay height obtained with the simulations increases logarithmically with *u*_{∗}. At low friction velocities, the high rate of increase in *δ*_{q} with *u*_{∗} suggests that the computed vertical ejection velocity scales with *u*_{∗}, as proposed by Owen (1964) and hypothesized by Nishimura and Hunt (2000) (see discussion in Sects. 2.3.1 and 2.3.3). In Fig. 5, we present the computed vertical profiles of the particle streamwise velocity (Fig. 5a) and of the vertical velocity of upward-moving particles (Fig. 5b). In the insets, the respective near-surface values are presented as a function of the friction velocity: ${\stackrel{\mathrm{\u203e}}{v}}_{x}^{\mathrm{e},\mathrm{i}}$ [m s^{−1}] and ${\stackrel{\mathrm{\u203e}}{v}}_{z}^{\mathrm{e}}$ [m s^{−1}] can be regarded as the particle slip velocity and the average vertical ejection velocity, respectively. For friction velocities ranging from 0.15 to 0.3 m s^{−1}, the particle streamwise velocity increases with height and with the friction velocity. The same trend is found for the vertical velocity of upward-moving particles in the first 2 cm above the surface. At these low values of *u*_{∗}, the average vertical ejection velocity and particle slip velocity also increase with the rise in the friction velocity. Therefore, the obtained increase in decay height as the friction velocity rises from 0.15 to 0.3 m s^{−1} is indeed accompanied by an increase in the vertical ejection velocity.

Conversely, a different trend is found for *u*_{∗} ranging between 0.39 and 0.7 m s^{−1}. For instance, an increase in the streamwise and vertical velocities with the increase in the friction velocity is only seen above a height of approximately 1.5 cm above the surface. Below this height, the particle velocity profile is mainly invariant with respect to *u*_{∗}, which is in agreement with several saltation models (e.g., Kok and Renno, 2009; Niiya and Nishimura, 2022) and measurements carried out over sand (Creyssels et al., 2009; Ho et al., 2011), discussed in Sects. 2.2.2 and 2.3.2. From the insets of Fig. 5, it can be seen that the near-surface velocity actually decreases with the increase in the friction velocity.

In the model, while the average ejection velocity of splashed particles scales with the particle impact velocity, the average ejection velocity of aerodynamically entrained particles scales with the surface friction velocity, ${u}_{\ast ,\mathrm{s}}$. The average vertical mass flux of aerodynamically entrained and splashed particles is presented in Fig. 6a, and the average surface friction velocity is presented in Fig. 6b, as a function of the friction velocity. Both quantities are computed by performing a surface average over the erodible bed and a time average over the last 100 s of the simulations. In Fig. 6a, it is seen that the vertical mass flux of aerodynamically entrained particles is greater than the vertical mass flux of particles ejected via splash for *u*_{∗} equal to 0.15 and 0.23 m s^{−1}. For *u*_{∗} = 0.3 m s^{−1}, splash entrainment is more significant than aerodynamic entrainment, but they are both of the same order of magnitude. For these three simulations, in which aerodynamic entrainment plays an important role during steady-state saltation, the average surface friction velocity does not deviate significantly from the friction velocity outside the saltation layer (Fig. 6b). Therefore, for *u*_{∗} ranging from 0.15 to 0.3 m s^{−1}, the increase in the vertical ejection velocity with the rise in *u*_{∗} is a direct consequence of the respective increase in the surface friction velocity.

For friction velocities greater than 0.3 m s^{−1}, the vertical mass flux of splashed particles is at least 1 order of magnitude greater than the vertical mass flux of aerodynamically entrained particles (Fig. 6a). In addition, the surface friction velocity deviates considerably from the friction velocity outside the saltation layer (Fig. 6b). It reaches a value smaller than *u*_{∗}, which decreases slightly with the rise in the friction velocity.

Even though the average ejection velocity of splashed particles scales with the velocity of the impacting particle and not directly with ${u}_{\ast ,\mathrm{s}}$, the computed particle impact velocity is highly correlated with the near-surface wind velocity (and, therefore, with the surface friction velocity) because of the exchange of momentum between the fluid and the particles close to the surface. The vertical profiles of the streamwise wind velocity are presented in Fig. 7. They are computed by averaging the streamwise wind velocity along horizontal planes and over the last 100 s of the simulations. Overall, similar features are seen in comparison to Fig. 5a: for friction velocities ranging from 0.15 to 0.3 m s^{−1}, for which aerodynamic entrainment plays an important role, the streamwise wind velocity increases with height and friction velocity; however, at greater values of *u*_{∗}, for which particle entrainment is dominated by splashed particles, the streamwise wind velocity profiles obtained during steady-state saltation exhibit a focus point at approximately 1 cm above the surface. As discussed in Sect. 2.3.2, the existence of a focus point was observed in wind tunnel experiments of sand saltation (Bagnold, 1941; Li and McKenna Neuman, 2012) and is supported by numerical models that represent the surface processes of rebound and splash (Kok et al., 2012; Durán et al., 2011). Due to the focus point, at high friction velocities, the near-surface wind velocity during steady-state saltation is approximately invariant with respect to *u*_{∗}. Nevertheless, below the focus point, the streamwise wind velocity decreases slightly as *u*_{∗} increases (see zoom-in to the near-surface region in the inset of Fig. 7). Below the first grid point above the surface (∼ 5 mm), the wind velocity is assumed to follow a logarithmic profile characterized by a constant roughness length. Therefore, the obtained increase in streamwise wind velocity at the first grid point above the surface as the friction velocity rises from 0.15 to 0.3 m s^{−1} implies an increase in the surface friction velocity with respect to *u*_{∗} (Fig. 6b). Conversely, the slight decrease in streamwise wind velocity at the first grid point above the surface as *u*_{∗} rises above 0.39 m s^{−1} justifies the slight decrease in ${u}_{\ast ,\mathrm{s}}$ as the friction velocity increases from 0.39 to 0.7 m s^{−1}, observed in Fig. 6b. The lower the wind velocity in the vicinity of the bed, the lower the ability of the flow to accelerate the particles close to the surface. Therefore, the obtained decrease in near-surface particle velocity as the friction velocity rises above 0.39 m s^{−1} (insets in Fig. 5) is due to the respective decrease in streamwise wind velocity below the focus point with respect to *u*_{∗}.

At friction velocities greater than 0.3 m s^{−1}, the current numerical model confirms neither the scaling of the ejection velocity with the friction velocity outside the saltation layer nor the scaling of the decay height with the vertical ejection velocity. In fact, even though the vertical ejection velocity decreases with the rise in the friction velocity for *u*_{∗} > 0.3 m s^{−1} (inset in Fig. 5b), the decay height increases monotonically with *u*_{∗} (Fig. 4). Therefore, the model proposed by Nishimura and Hunt (2000) might not be appropriate to parameterize the decay height during splash-dominated saltation.

The obtained evolution of ${u}_{\ast ,\mathrm{s}}$ at low friction velocities presented in Fig. 6b contrasts with the results of sand saltation models. For instance, Kok and Renno (2009) neglected the surface process of aerodynamic entrainment and found ${u}_{\ast ,\mathrm{s}}$ to decrease monotonically with the increase in the friction velocity. Sand particles have a higher density than snow particles, which implies a higher fluid threshold. In addition, as discussed in Sect. 2.2.1, the surface shear stress is expected to be lower than the fluid threshold, which restricts the occurrence of aerodynamic entrainment during steady-state saltation, even at low friction velocities.

In contrast, similar features are seen in some snow saltation models in which aerodynamic entrainment is represented. For example, Nemoto and Nishimura (2004) found the surface friction velocity to increase monotonically with the rise in the friction velocity from 0.23 to 0.39 m s^{−1}. The same increase in ${u}_{\ast ,\mathrm{s}}$ is seen in Fig. 6b for this range of friction velocities. The values of ${u}_{\ast ,\mathrm{s}}$ found by Nemoto and Nishimura (2004) are greater than the specified fluid threshold friction velocity, which suggests that aerodynamic entrainment also plays a role in their steady-state conditions. The process of aerodynamic entrainment is also taken into account in the numerical model of Doorschot and Lehning (2002). In their model, the parameterizations for aerodynamic entrainment and rebound are similar to those described in Sect. 3.1, but the saltation system is – for simplicity – assumed to be either entirely composed of aerodynamically entrained particles or of particles that continuously rebound from the surface (in their model, no distinction is made between rebound and splash). Similarly to the results presented in Fig. 6, Doorschot and Lehning (2002) found saltation to be dominated by aerodynamic entrainment at low friction velocities and to be dominated by rebounding and/or splashed particles at high friction velocities. In addition, at low friction velocities, the surface friction velocity increases with the rise in *u*_{∗}, and a negligible deviation is seen between the values of ${u}_{\ast ,\mathrm{s}}$ and *u*_{∗} due to the low number of particles aloft. At high friction velocities, ${u}_{\ast ,\mathrm{s}}$ decreases slightly with the rise in *u*_{∗}. In their model, the friction velocity at which the transition from aerodynamic entrainment to rebound and/or splash occurs depends primarily on the assumed restitution coefficient. This coefficient is expected to be highly dependent on the snow surface characteristics, which highlights the importance of the snow cover properties in the correct description of snow saltation and its scaling laws.

In the snow saltation model of Niiya and Nishimura (2022), ${u}_{\ast ,\mathrm{s}}$ is found to decrease slightly with the increase in the friction velocity from 0.24 to 0.3 m s^{−1}. This differs from the increase in ${u}_{\ast ,\mathrm{s}}$ found in Fig. 6b for the same range of friction velocities. In their work, the fluid threshold friction velocity is assumed equal to 0.24 m s^{−1}, and the surface friction velocity is found to be approximately equal to 0.2 m s^{−1}. In comparison to the current simulations, Niiya and Nishimura (2022) assumed a higher value for the fluid threshold friction velocity (0.24 versus 0.154 m s^{−1}) and a lower value for the roughness length (10^{−5} versus 10^{−4} m). When decreasing *z*_{o} from 10^{−4} to 10^{−5} m in the LES-LSM, we obtain a significant decrease in the surface friction velocity during splash-dominated saltation. In this way, a decrease in *z*_{o} and an increase in *τ*_{ft} reduce the range of friction velocities for which saltation is controlled by aerodynamic entrainment, which might justify the seemingly different results obtained with both models. In addition to the roughness length and the fluid threshold, other parameters that characterize the snow surface and highly influence the surface processes of aerodynamic entrainment, rebound and splash are expected to influence the importance of aerodynamic entrainment during steady-state saltation. The narrower the range of friction velocities for which aerodynamic entrainment plays an important role, the narrower the range of friction velocities for which the decay height is well approximated by the expression of Nishimura and Hunt (2000).

The modeled saltation system is further investigated by the analysis of the particle hop height and length. The average hop height and length of saltating particles are computed from the analysis of all particle hops during the last 1 s of simulations. During this time period, the location of all particles aloft was outputted by the model at very high frequency. In this analysis, a particle hop is defined as a fraction of a particle trajectory limited by two consecutive local minimums located close to the surface. The hop height is given by the maximum height attained by the particle along the respective hop. An example of a particle trajectory comprising several hops is presented in Fig. 8a. The particles seldom reach the surface because the surface processes of rebound and splash are assumed to take place as soon as the particles reach a height lower than $z=\mathrm{4}\stackrel{\mathrm{\u203e}}{d}$ = 0.14 cm.

The resultant average hop height and length are presented in Fig. 8b as a functions of the friction velocity. The average hop height increases with the rise in the friction velocity until *u*_{∗} = 0.5 m s^{−1}. It increases at a higher rate for *u*_{∗} ranging from 0.15 to 0.3 m s^{−1}, which is in agreement with the evolution found for the vertical ejection velocity and the decay height (Figs. 4 and 5b). As *u*_{∗} increases above 0.3 m s^{−1} and saltation becomes splash-dominated, the average hop height increases at a lower rate with the rise in the friction velocity. It reaches a plateau at *u*_{∗} = 0.5 m s^{−1} of approximately 2.8 cm. For *u*_{∗} > 0.3 m s^{−1}, even though the vertical ejection velocity decreases with the rise in the friction velocity (Fig. 5b), 〈*h*_{p}〉 is approximately invariant with respect to *u*_{∗}. This apparent discrepancy might be related to the fact that the hop height does not depend solely on the vertical ejection velocity but also on the vertical component of the aerodynamic force applied to a particle along its trajectory.

As discussed in Sects. 2.3.1 and 2.3.3, the decay height is expected to scale with the average hop height of particles in saltation (Nishimura and Hunt, 2000; Gordon et al., 2009; Martin and Kok, 2017). However, in the model, a linear relationship between 〈*h*_{p}〉 and *δ*_{q} is not found for *u*_{∗} > 0.3 m s^{−1}. The modeled saltation system is not only composed of particles undergoing ballistic trajectories, but it also includes particles traveling without regular contact with the surface in the vicinity of the saltation–suspension interface. In this context, it is reasonable to assume that *δ*_{q} does not depend entirely on the hop height of ballistic hops but rather on the particle motion of all particles that contribute to the particle mass flux in the saltation layer.

The average particle hop length, 〈*l*_{p}〉, increases continuously with the rise in the friction velocity. Similarly to the average hop height, it increases at a higher rate when *u*_{∗} varies between 0.15 and 0.3 m s^{−1} than when *u*_{∗} is greater than 0.3 m s^{−1} (Fig. 8b). In this way, when saltation is dominated by splash, the particle hop length can be approximated as invariant with respect to *u*_{∗}, which is in agreement with the saltation model of Ungar and Haff (1987) discussed in Sect. 2.2.3. If we define the impact threshold friction velocity as the minimum surface friction velocity for which saltation is dominated by splash (in this case, ${u}_{\ast ,\text{it}}$ ≈ 0.31 m s^{−1}), the studied values of friction velocity are within the low-velocity regime defined by Durán et al. (2011) when modeling sand saltation (${u}_{\ast}/{u}_{\ast ,\text{it}}<\mathrm{4}$). At these friction velocities, Durán et al. (2011) expect saltation to be dominated by splash and the average hop length to be invariant with respect to *u*_{∗}. Indeed, when saltation is dominated by splash, the current results are in agreement with those from Durán et al. (2011). This more restrictive definition of the impact threshold friction velocity (recall the definition of ${u}_{\ast ,\text{it}}$ presented in Sect. 2.2.1) is considered in some saltation models (e.g., Comola et al., 2022).

If 〈*l*_{p}〉 is mainly invariant with respect to *u*_{∗}, the transport rate is expected to scale with ${u}_{\ast}^{\mathrm{2}}$. The transport rate is presented in Fig. 9 as a function of the friction velocity. It is computed by integrating the particle mass flux along the height, from the surface to a height of 15 cm. The fits between the computed values obtained for *u*_{∗} > 0.3 m s^{−1} and a quadratic function, as well as a cubic function, are also presented. Indeed, the transport rate obtained during splash-dominated saltation is in good agreement with a quadratic curve, but the root-mean-square error (RMSE) of the cubic function is only slightly higher. When considering the full range of friction velocities, a better agreement is actually found with the cubic function. These results highlight the need to acquire mass flux measurements for a wide range of friction velocities and snow surface characteristics to fully evaluate the scaling of the transport rate with the friction velocity.

The two saltation regimes obtained with the model (dominated by aerodynamic entrainment or splash) are illustrated in Fig. 10. In summary, as the friction velocity increases, the saltation system evolves from an aerodynamic-entrainment-dominated system to one dominated by splashed particles. In the first regime, the particle streamwise velocity increases with the rise in the friction velocity at all heights, and the average hop height of saltating particles (represented by the line enclosing the ballistic trajectories) increases as well. As the friction velocity reaches higher values, the average hop height stops increasing significantly with the rise in the friction velocity. In fact, once saltation is dominated by splash, the average hop height is approximately invariant with respect to *u*_{∗}. In addition, the particle streamwise velocity is also mainly invariant with respect to *u*_{∗} in the near-surface region where most saltating particles are. Considering the alternative definition for the impact threshold friction velocity previously mentioned, one can state that these two saltation regimes are typical of saltation systems characterized by an impact threshold friction velocity, ${u}_{\ast ,\text{it}}$, greater than the fluid threshold friction velocity, ${u}_{\ast ,\text{ft}}$. As previously discussed, this is not characteristic of sand saltation, but it is probably representative of snow saltation for some snow surface characteristics. In Fig. 10, the line indicating the average hop height of saltating particles separates the regions where particle motion is controlled by wind–particle interaction only (above) and where it is highly dependent on the particle–bed interaction (below).

The conceptual model of snow saltation illustrated in Fig. 10 neglects the effect of large turbulent eddies that exist over snow-covered regions subjected to high friction velocities. Indeed, these large eddies will enhance the entrainment of particles in suspension and perhaps limit or suppress the motion of particles in saltation. This limiting case goes beyond the scope of this study and needs further investigation. It can neither be evaluated with the computational domain considered in this analysis nor with wind tunnel experiments. In addition, it represents a challenge for the instruments currently available to assess snow saltation dynamics.

The importance of drifting snow events to the mass and energy balances of snow-covered regions has long been acknowledged by the scientific community. Several atmospheric and snow models are currently enriched with drifting snow schemes that represent the eolian transport of snow and quantify the induced changes in snow height and the amount of snow sublimation. Nevertheless, the correct prediction of these quantities is still a challenge. This is due to the complexity of the phenomenon and the uncertainties in the parameters that control drifting snow transport and sublimation but also due to inaccuracies in the drifting snow schemes themselves.

In this work, we focus on an important part of drifting snow schemes that is typically overlooked: the saltation model. By comparing the saltation models used in the drifting snow schemes of RACMO, MAR, Meso-NH and CRYOWRF (Amory et al., 2021; Vionnet et al., 2014; van Wessem et al., 2018; Sharma et al., 2023) with the different theoretical, experimental and numerical studies on sand and snow saltation, we conclude that the parameterizations employed are not fully aligned with the current understanding of snow saltation dynamics. In addition, by performing numerical simulations with an LES-based model, we verify the consistency of some parameterizations and show the relationship between them and the prevalence of aerodynamically entrained or splashed particles.

In the saltation model of RACMO and MAR, the particle mass flux in the saltation layer is assumed uniform in height (Pomeroy and Gray, 1990). This contrasts with the well-documented exponential decay found in wind tunnel and field observations of sand and snow saltation, as well as in numerical simulations. In addition, the calculation of the transport rate and of the particle concentration at the top of the saltation layer is based on the assumption that the saltation layer height is proportional to ${u}_{\ast}^{\mathrm{2}}$ and that the average particle streamwise velocity in the saltation layer is invariant with respect to *u*_{∗}. Even though both scaling laws can be obtained from theoretical, numerical or experimental studies, depending on the friction velocity and on the snow surface characteristics, these two assumptions are contradictory. While the first is based on the idea that the near-surface particle velocity scales with *u*_{∗}, the latter is based on the notion that the near-surface particle velocity is invariant with respect to *u*_{∗}. As a result, the transport rate is found to scale linearly with *u*_{∗}, which contrasts with most saltation models and measurements available in the literature (e.g., Doorschot and Lehning, 2002; Sørensen, 2004; Melo et al., 2022). Conversely, the saltation models used in Meso-NH and CRYOWRF assume that the particle mass flux profile follows an exponential decay. However, when computing the decay height and the particle concentration, which are a function of the height of the saltation layer and of the average particle streamwise velocity, they consider the same contradictory parameterizations regarding the near-surface particle velocity.

The available experimental and numerical studies show that the near-surface particle velocity is invariant with respect to *u*_{∗} when saltation develops over a sand-covered surface. In this case, the average hop height of saltating particles, the height of the saltation layer, and the decay height of the exponential mass flux and concentration profiles are found to be equally invariant with respect to *u*_{∗}. However, when saltation develops over snow-covered surfaces, the scaling of these different quantities with *u*_{∗} is less clear. While early wind tunnel studies reveal that the decay height of the mass flux profile scales with ${u}_{\ast}^{\mathrm{2}}$ (and, therefore, that the near-surface particle velocity scales with *u*_{∗}) (Sugiura et al., 1998), recent experiments do not reveal a clear scaling of the decay height with the friction velocity, and they show that the scaling of the near-surface particle velocity with *u*_{∗} depends on the snow surface characteristics (Aksamit and Pomeroy, 2016; Gordon et al., 2009).

From numerical simulations, we show that the scaling of the near-surface particle velocity, average hop height and decay height are a function of the saltation regime: when saltation is dominated by aerodynamic entrainment, the near-surface particle velocity and the average hop height increase with the friction velocity; when saltation is dominated by splash, these two quantities are mainly invariant with respect to *u*_{∗}. For the full range of friction velocities studied, the decay height is found to increase logarithmically with *u*_{∗}. Therefore, the rate of increase in the decay height with *u*_{∗} decreases as *u*_{∗} rises. Thus, when saltation is dominated by splash, the decay height increases at a much lower rate than what is obtained when saltation is dominated by aerodynamic entrainment. When saltation is dominated by splash, the decay height significantly deviates from the evolution proposed by Nishimura and Hunt (2000), which is considered in the drifting snow schemes of Meso-NH and CRYOWRF. The effect of the saltation regime on the saltation dynamics and its scaling laws might partly justify the different results obtained by Sugiura et al. (1998) in a wind tunnel in comparison to the results obtained by Gordon et al. (2009) and Aksamit and Pomeroy (2016) in the field. Indeed, the limited size of a wind tunnel might restrict the development of a fully developed saltation system dominated by splash. This increases the relative importance of aerodynamic entrainment and might justify the observed scaling of the decay height with ${u}_{\ast}^{\mathrm{2}}$ obtained by Sugiura et al. (1998). Other aspects, such as the height at which the particle concentration and mass flux are measured and the snow surface characteristics, are also expected to influence the decay height and its scaling with *u*_{∗}. Further investigation is needed to fully understand the impact of these different aspects on the decay height.

In general, the saltation models used in RACMO, MAR, Meso-NH and CRYOWRF have a poor representation of the effect of the snow surface characteristics on the different quantities of interest. The snow surface characteristics have a direct effect on one single parameter – the fluid threshold – which is used as an approximation for the surface shear stress. Therefore, they influence the transport rate and the average particle streamwise velocity. However, even though both the fluid threshold and the surface shear stress are found to increase as the grain diameter increases, the opposite is found for the particle streamwise velocity from previous works. In addition, the numerical simulations suggest that the surface shear stress might not be well approximated by the fluid threshold and that it increases with the rise in the friction velocity when saltation is dominated by aerodynamic entrainment. Therefore, improvements are needed to correctly take into account the effect of the snow surface characteristics on the snow saltation parameterizations. For instance, the coefficients in the expression of Sørensen (2004) used in Meso-NH and CRYOWRF to compute the transport rate are assumed constant. However, they are expected to vary with the snow surface characteristics. Extensive measurement campaigns are needed in order to fully understand and quantify the effect of the snow surface on the saltation dynamics. In addition, it is essential that future experimental studies describe the characteristics of the snow surface, the wind field and the particles in saltation.

In this work, we have mainly analyzed the parameterizations used to compute the transport rate and the vertical profiles of particle mass flux, concentration and streamwise velocity. However, the description of snow saltation in a drifting snow scheme is also dependent on the parameterizations used for the fluid threshold, the particle size distribution and the implementation of the lower boundary condition for the snow suspension scheme. Therefore, a complete assessment of a saltation model must take these parameterizations into account. In addition, the height at which the lower boundary condition is specified must be in agreement with the assumed parameterizations. The actual value used in drifting snow schemes (which is typically given by the height of the saltation layer) is rather irrelevant. However, one must guarantee that all parameterizations that constitute the saltation model are valid at that height. For instance, when saltation is dominated by splash, the presented simulation results show that the average particle streamwise velocity is only invariant with respect to *u*_{∗} in the first 1–2 cm above the surface. Even though saltating particles can reach higher elevations, the use of the previous assumption supports the application of the lower boundary condition close to the surface. Indeed, all drifting snow schemes studied set the lower boundary condition at a height relatively close to the surface: typically, its value increases with *u*_{∗}, reaching a value of 5–6 cm for *u*_{∗} = 0.8 m s^{−1} (Pomeroy and Gray, 1990; Pomeroy and Male, 1992). However, the particle size distribution at the lower boundary is generally assumed independent of the snow surface characteristics (Vionnet et al., 2014; Sharma et al., 2023), which is not a reasonable approximation close to the surface but only at higher elevations, close to the saltation–suspension interface (8–10 cm). The lower boundary condition of snow suspension schemes can be specified at this location. However, the parameterization for the particle streamwise velocity would have to be adjusted, in particular its scaling with *u*_{∗}. Even though this would add additional unknown parameters to the model, the use of an advanced parameterization for the particle streamwise velocity would allow the description of saltation systems dominated by aerodynamic entrainment. Nevertheless, when setting the lower boundary at a higher elevation, one must bear in mind that the particle mass flux profile is expected to deviate from an exponential decay at the saltation–suspension interface (e.g., Nishimura and Nemoto, 2005). Taking into account the advantages and disadvantages, this choice is left to the modeler.

The snow saltation parameterizations studied in this work are based on the steady-state saltation assumption. Even though this approximation is widely used, it does not fully represent the phenomenon of snow saltation in natural environments. Non-stationary wind conditions (Aksamit and Pomeroy, 2018) and intermittent saltation (Paterna et al., 2016; Comola et al., 2019c) have been recently studied, but further investigation is needed to incorporate these findings in simple parameterizations.

The correct prediction of snow transport and sublimation at large scales has a variety of challenges. One of them is the representation of snow saltation by means of simple parameterizations. Snow saltation is a complex phenomenon, challenging to measure and model at both small and large scales. The full understanding of the processes involved is a multi-disciplinary endeavor, which requires contributions from both atmospheric and snow scientists. By means of a thorough and careful analysis of previous developments and new simulations, this work offers a new insight into the saltation system and presents some pathways to improve the representation of snow saltation in atmospheric models.

The code version used in the simulations, the results presented, and the post-processing scripts used to analyze the simulation outputs and produce the figures are available at the EnviDat repository (https://doi.org/10.16904/envidat.479, Melo et al., 2024). The code is also available at the institutional GitLab repository (https://gitlabext.wsl.ch/atmospheric-models/les-lsm (last access: 30 January 2024) – see tag “r2024”).

DBM performed the simulations, implemented the post-processing scripts and wrote the literature review. DBM, AS and ML contributed to the analysis of the results. AS and DBM contributed to the code development and its maintenance. ML supervised the work and secured its funding. DBM prepared the article with contributions from all co-authors.

The contact author has declared that none of the authors has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

The authors acknowledge the Swiss National Supercomputing Centre (CSCS) for providing the computational resources (project nos. s1031 and s1115). Pedro Cabral is acknowledged for the numerous in-depth discussions about physics and fluid dynamics (which contributed to the presented literature review and subsequent analysis), for preparing Figs. 1 and 10, and for the detailed revision of the text. Varun Sharma is acknowledged for his support in the early stages of this work and for several discussions that have inspired part of the work developed. Raleigh L. Martin, Kouichi Nishimura, Nikolas O. Aksamit, Masaki Nemoto, and Konosuke Sugiura are acknowledged for promptly providing extra clarifications about their works. Finally, the authors thank the two reviewers for their constructive comments and suggestions, which have improved the clarity and completeness of this article.

This research has been supported by the Swiss National Science Foundation (SNSF) (grant no. 179130).

This paper was edited by Masashi Niwano and reviewed by Nikolas Aksamit and one anonymous referee.

Aksamit, N. O. and Pomeroy, J. W.: Near-surface snow particle dynamics from particle tracking velocimetry and turbulence measurements during alpine blowing snow storms, The Cryosphere, 10, 3043–3062, https://doi.org/10.5194/tc-10-3043-2016, 2016. a, b, c, d, e, f, g, h, i, j

Aksamit, N. O. and Pomeroy, J. W.: Scale Interactions in Turbulence for Mountain Blowing Snow, J. Hydrometeorol., 19, 305–320, https://doi.org/10.1175/JHM-D-17-0179.1, 2018. a, b

Albertson, J. D. and Parlange, M. B.: Natural integration of scalar fluxes from complex terrain, Adv. Water Resour., 23, 239–252, https://doi.org/10.1016/S0309-1708(99)00011-1, 1999. a

Amory, C., Trouvilliez, A., Gallée, H., Favier, V., Naaim-Bouvet, F., Genthon, C., Agosta, C., Piard, L., and Bellot, H.: Comparison between observed and simulated aeolian snow mass fluxes in Adélie Land, East Antarctica, The Cryosphere, 9, 1373–1383, https://doi.org/10.5194/tc-9-1373-2015, 2015. a

Amory, C., Kittel, C., Le Toumelin, L., Agosta, C., Delhasse, A., Favier, V., and Fettweis, X.: Performance of MAR (v3.11) in simulating the drifting-snow climate and surface mass balance of Adélie Land, East Antarctica, Geosci. Model Dev., 14, 3487–3510, https://doi.org/10.5194/gmd-14-3487-2021, 2021. a, b, c, d, e, f, g, h, i

Anderson, R. S. and Haff, P. K.: Wind modification and bed response during saltation of sand in air, in: Aeolian Grain Transport 1, edited by: Barndorff-Nielsen, O. E. and Willetts, B. B., Springer Vienna–New York, pp. 21–51, ISBN 978-3-211-82269-2, https://doi.org/10.1007/978-3-7091-6706-9, 1991. a, b

Anderson, R. S., Sørensen, M., and Willetts, B. B.: A review of recent progress in our understanding of aeolian sediment transport, in: Aeolian Grain Transport 1, edited by: Barndorff-Nielsen, O. E. and Willetts, B. B., Springer Vienna–New York, pp. 1–19, ISBN 978-3-211-82269-2, https://doi.org/10.1007/978-3-7091-6706-9, 1991. a

Araoka, K. and Maeno, N.: Dynamical behaviors of snow particles in the saltation layer, in: Proceedings of the Third Symposium on Polar Meteorology and Glaciology, 13–14 January 1981, National Institute of Polar Research, Tokyo, vol. 19, pp. 253–263, 1981. a

Bagnold, R. A.: The Physics of Blown Sand and Desert Dunes, Dover Publications, Mineola, New York, ISBN-13 978-0-486-43931-0, 1941. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o

Bagnold, R. A.: The nature of saltation and of '1bed-load' transport in water, P. R. Soc. Lond. A, 332, 473–504, https://doi.org/10.1098/rspa.1973.0038, 1973. a, b, c, d

Bartelt, P. and Lehning, M.: A physical SNOWPACK model for the Swiss avalanche warning: Part I: numerical model, Cold Reg. Sci. Technol., 35, 123–145, https://doi.org/10.1016/S0165-232X(02)00074-5, 2002. a

Bauer, B. O. and Davidson-Arnott, R. G. D.: Aeolian particle flux profiles and transport unsteadiness, J. Geophys. Res.-Earth, 119, 1542–1563, https://doi.org/10.1002/2014JF003128, 2014. a, b

Bernhardt, M., Zängl, G., Liston, G. E., Strasser, U., and Mauser, W.: Using wind fields from a high-resolution atmospheric model for simulating snow dynamics in mountainous terrain, Hydrol. Process., 23, 1064–1075, https://doi.org/10.1002/hyp.7208, 2009. a

Bintanja, R.: The contribution of snowdrift sublimation to the surface mass balance of Antarctica, Ann. Glaciol., 27, 251–259, https://doi.org/10.3189/1998AoG27-1-251-259, 1998. a

Bou-Zeid, E., Meneveau, C., and Parlange, M.: A scale-dependent Lagrangian dynamic model for large eddy simulation of complex turbulent flows, Phys. Fluids, 17, 025105, https://doi.org/10.1063/1.1839152, 2005. a

Budd, W. F., Dingle, W. R. J., and Radok, U.: The Byrd Snow Drift Project: Outline and Basic Results, in: Studies in Antarctic Meteorology, vol. 9 of Antartic Research Series, edited by: Rubin, M. J., pp. 71–134, https://doi.org/10.1029/AR009p0071, 1966. a

Clifton, A., Rüedi, J.-D., and Lehning, M.: Snow saltation threshold measurements in a drifting-snow wind tunnel, J. Glaciol., 52, 585–596, https://doi.org/10.3189/172756506781828430, 2006. a

Comola, F.: Stochastic modeling of snow transport and hydrologic response in alpine terrain, PhD thesis, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland, https://doi.org/10.5075/epfl-thesis-7532, 2017. a

Comola, F. and Lehning, M.: Energy- and momentum-conserving model of splash entrainment in sand and snow saltation, Geophys. Res. Lett., 44, 1601–1609, https://doi.org/10.1002/2016GL071822, 2017. a, b, c, d

Comola, F., Kok, J. F., Gaume, J., Paterna, E., and Lehning, M.: Fragmentation of wind-blown snow crystals, Geophys. Res. Lett., 44, 4195–4203, https://doi.org/10.1002/2017GL073039, 2017. a

Comola, F., Gaume, J., Kok, J. F., and Lehning, M.: Cohesion-Induced Enhancement of Aeolian Saltation, Geophys. Res. Lett., 46, 5566–5574, https://doi.org/10.1029/2019GL082195, 2019a. a

Comola, F., Giometto, M. G., Salesky, S. T., Parlange, M. B., and Lehning, M.: Preferential Deposition of Snow and Dust Over Hills: Governing Processes and Relevant Scales, J. Geophys. Res.-Atmos., 124, 7951–7974, https://doi.org/10.1029/2018JD029614, 2019b. a, b, c

Comola, F., Kok, J. F., Chamecki, M., and Martin, R. L.: The Intermittency of Wind-Driven Sand Transport, Geophys. Res. Lett., 46, 13430–13440, https://doi.org/10.1029/2019GL085739, 2019c. a

Comola, F., Kok, J. F., Lora, J. M., Cohanim, K., Yu, X., He, C., McGuiggan, P., Hörst, S. M., and Turney, F.: Titan's Prevailing Circulation Might Drive Highly Intermittent, Yet Significant Sediment Transport, Geophys. Res. Lett., 49, e2022GL097913, https://doi.org/10.1029/2022GL097913, 2022. a, b

Creyssels, M., Dupont, P., Moctar, A. O. E., Valance, A., Cantat, I., Jenkins, J. T., Pasini, J. M., and Rasmussen, K. R.: Saltating particles in a turbulent boundary layer: experiment and theory, J. Fluid Mech., 625, 47–74, https://doi.org/10.1017/S0022112008005491, 2009. a, b, c, d, e, f, g, h, i, j, k, l

Dai, X. and Huang, N.: Numerical simulation of drifting snow sublimation in the saltation layer, Sci. Rep.-UK, 4, 6611, https://doi.org/10.1038/srep06611, 2014. a

Das, R. K., Datt, P., and Acharya, A.: An assessment of the FlowCapt acoustic sensor for measuring snowdrift in the Indian Himalayas, J. Earth Syst. Sci., 121, 1483–1491, https://doi.org/10.1007/s12040-012-0234-2, 2012. a

Déry, S. J., Taylor, P., and Xiao, J.: The Thermodynamic Effects of Sublimating, Blowing Snow in the Atmospheric Boundary Layer, Bound.-Lay. Meteorol., 89, 251–283, https://doi.org/10.1023/A:1001712111718, 1998. a, b

Doorschot, J. J. J. and Lehning, M.: Equilibrium Saltation: Mass Fluxes, Aerodynamic Entrainment, and Dependence on Grain Properties, Bound.-Lay. Meteorol., 104, 111–130, https://doi.org/10.1023/A:1015516420286, 2002. a, b, c, d, e, f, g

Durán, O., Claudin, P., and Andreotti, B.: On aeolian transport: Grain-scale interactions, dynamical mechanisms and scaling laws, Aeolian Res., 3, 243–270, https://doi.org/10.1016/j.aeolia.2011.07.006, 2011. a, b, c, d, e, f, g, h, i, j, k

Dyunin, A. K.: Fundamentals of the Mechanics of Snow Storms, Physics of Snow and Ice: proceedings, 1, in: Physics of Snow and Ice: proceedings, 1, 1065–1073, edited by: Oura, H., International Conference on Low Temperature Science, Sapporo, 14–19 August 1966, Hokkaido University, Institute of Low Temperature Science, 1967. a

Dyunin, A. K. and Kotlyakov, V. M.: Redistribution of snow in the mountains under the effect of heavy snow-storms, Cold Reg. Sci. Technol., 3, 287–294, https://doi.org/10.1016/0165-232X(80)90035-X, 1980. a

Filhol, S. and Sturm, M.: Snow bedforms: A review, new data, and a formation model, J. Geophys. Res.-Earth, 120, 1645–1669, https://doi.org/10.1002/2015JF003529, 2015. a

Gallée, H.: Simulation of blowing snow over the antarctic ice sheet, Ann. Glaciol., 26, 203–206, https://doi.org/10.3189/1998AoG26-1-203-206, 1998. a

Gallée, H., Guyomarc'h, G., and Brun, E.: Impact Of Snow Drift On The Antarctic Ice Sheet Surface Mass Balance: Possible Sensitivity To Snow-Surface Properties, Bound.-Lay. Meteorol., 99, 1–19, https://doi.org/10.1023/A:1018776422809, 2001. a, b

Gauer, P.: Blowing and drifting snow in Alpine terrain: numerical simulation and related field measurements, Ann. Glaciol., 26, 174–178, https://doi.org/10.3189/1998AoG26-1-174-178, 1998. a

Gerber, F., Sharma, V., and Lehning, M.: CRYOWRF – Model Evaluation and the Effect of Blowing Snow on the Antarctic Surface Mass Balance, J. Geophys. Res.-Atmos., 128, e2022JD037744, https://doi.org/10.1029/2022JD037744, 2023. a, b

Giometto, M. G., Christen, A., Meneveau, C., Fang, J., Krafczyk, M., and Parlange, M. B.: Spatial Characteristics of Roughness Sublayer Mean Flow and Turbulence Over a Realistic Urban Surface, Bound.-Lay. Meteorol., 160, 425–452, https://doi.org/10.1007/s10546-016-0157-6, 2016. a

Giometto, M. G., Christen, A., Egli, P. E., Schmid, M. F., Tooke, R. T., Coops, N. C., and Parlange, M. B.: Effects of trees on mean wind, turbulence and momentum exchange within and above a real urban environment, Adv. Water Resour., 106, 154–168, https://doi.org/10.1016/j.advwatres.2017.06.018, 2017. a

Gordon, M., Savelyev, S., and Taylor, P. A.: Measurements of blowing snow, part II: Mass and number density profiles and saltation height at Franklin Bay, NWT, Canada, Cold Reg. Sci. Technol., 55, 75–85, https://doi.org/10.1016/j.coldregions.2008.07.001, 2009. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o

Groot Zwaaftink, C. D., Löwe, H., Mott, R., Bavay, M., and Lehning, M.: Drifting snow sublimation: A high-resolution 3-D model with temperature and moisture feedbacks, J. Geophys. Res.-Atmos., 116, D16107, https://doi.org/10.1029/2011JD015754, 2011. a, b

Groot Zwaaftink, C. D., Diebold, M., Horender, S., Overney, J., Lieberherr, G., Parlange, M. B., and Lehning, M.: Modelling Small-Scale Drifting Snow with a Lagrangian Stochastic Model Based on Large-Eddy Simulations, Bound.-Lay. Meteorol., 153, 117–139, https://doi.org/10.1007/s10546-014-9934-2, 2014. a

Guyomarc'h, G. and Mérindol, L.: Validation of an application for forecasting blowing snow, Ann. Glaciol., 26, 138–143, https://doi.org/10.3189/1998AoG26-1-138-143, 1998. a

Hancock, H., Eckerstorfer, M., Prokop, A., and Hendrikx, J.: Quantifying seasonal cornice dynamics using a terrestrial laser scanner in Svalbard, Norway, Nat. Hazards Earth Syst. Sci., 20, 603–623, https://doi.org/10.5194/nhess-20-603-2020, 2020. a

Ho, T. D., Valance, A., Dupont, P., and Ould El Moctar, A.: Scaling Laws in Aeolian Sand Transport, Phys. Rev. Lett., 106, 094501, https://doi.org/10.1103/PhysRevLett.106.094501, 2011. a, b, c, d, e, f, g, h, i

Jiménez, J.: Turbulent Flows over Rough Walls, Annu. Rev. Fluid Mech., 36, 173–196, https://doi.org/10.1146/annurev.fluid.36.050802.122103, 2004. a

Kawamura, R.: Study of sand movement by wind, Reports of Physical Sciences Research Institute of Tokyo University, NASA Technical Translation (NASA-TT-F-14215), National Aeronautics and Space Administration, Washington, D.C., 1972, 1951. a, b, c, d

Kobayashi, D.: Studies of Snow Transport In Low-Level Drifting Snow, Contributions from the Institute of Low Temperature Science, Institute of Low Temperature Science, Hokkaido University, A24, 1–58, 1972. a

Kok, J. F. and Renno, N. O.: A comprehensive numerical model of steady state saltation (COMSALT), J. Geophys. Res.-Atmos., 114, D17204, https://doi.org/10.1029/2009JD011702, 2009. a, b, c, d, e, f, g, h

Kok, J. F., Parteli, E. J. R., Michaels, T. I., and Karam, D. B.: The physics of wind-blown sand and dust, Rep. Prog. Phys., 75, 106901, https://doi.org/10.1088/0034-4885/75/10/106901, 2012. a, b, c, d, e, f, g, h

Lehning, M., Doorschot, J., and Bartelt, P.: A snowdrift index based on SNOWPACK model calculations, Ann. Glaciol., 31, 382–386, https://doi.org/10.3189/172756400781819770, 2000. a, b

Lehning, M., Löwe, H., Ryser, M., and Raderschall, N.: Inhomogeneous precipitation distribution and snow transport in steep terrain, Water Resour. Res., 44, W07404, https://doi.org/10.1029/2007WR006545, 2008. a, b

Lenaerts, J. T. M. and van den Broeke, M. R.: Modeling drifting snow in Antarctica with a regional climate model: 2. Results, J. Geophys. Res.-Atmos., 117, D05109, https://doi.org/10.1029/2010JD015419, 2012. a

Lenaerts, J. T. M., van den Broeke, M. R., Déry, S. J., van Meijgaard, E., van de Berg, W. J., Palm, S. P., and Sanz Rodrigo, J.: Modeling drifting snow in Antarctica with a regional climate model: 1. Methods and model evaluation, J. Geophys. Res.-Atmos., 117, D05108, https://doi.org/10.1029/2011JD016145, 2012. a, b, c, d, e, f, g, h

Li, B. and McKenna Neuman, C.: Boundary-layer turbulence characteristics during aeolian saltation, Geophys. Res. Lett., 39, L11402, https://doi.org/10.1029/2012GL052234, 2012. a, b, c

Li, L. and Pomeroy, J. W.: Estimates of Threshold Wind Speeds for Snow Transport Using Meteorological Data, J. Appl. Meteorol. Clim., 36, 205–213, https://doi.org/10.1175/1520-0450(1997)036<0205:EOTWSF>2.0.CO;2, 1997. a

Liston, G. and Sturm, M.: The role of winter sublimation in the Arctic moisture budget, Hydrol. Res., 35, 325–334, https://doi.org/10.2166/nh.2004.0024, 2004. a

Liston, G. E. and Sturm, M.: Winter Precipitation Patterns in Arctic Alaska Determined from a Blowing-Snow Model and Snow-Depth Observations, J. Hydrometeorol., 3, 646–659, https://doi.org/10.1175/1525-7541(2002)003<0646:WPPIAA>2.0.CO;2, 2002. a

Martin, R. L. and Kok, J. F.: Wind-invariant saltation heights imply linear scaling of aeolian saltation flux with shear stress, Science Advances, 3, e1602569, https://doi.org/10.1126/sciadv.1602569, 2017. a, b, c, d, e, f, g, h, i, j

Martin, R. L., Kok, J. F., Hugenholtz, C. H., Barchyn, T. E., Chamecki, M., and Ellis, J. T.: High-frequency measurements of aeolian saltation flux: Field-based methodology and applications, Aeolian Res., 30, 97–114, https://doi.org/10.1016/j.aeolia.2017.12.003, 2018. a

Melo, D. B., Sharma, V., Comola, F., Sigmund, A., and Lehning, M.: Modeling Snow Saltation: The Effect of Grain Size and Interparticle Cohesion, J. Geophys. Res.-Atmos., 127, e2021JD035260, https://doi.org/10.1029/2021JD035260, 2022. a, b, c, d, e, f, g, h, i, j, k, l

Melo, D. B., Sigmund, A., and Lehning, M.: Understanding snow saltation parameterizations – Reproducibility, EnviDat [code and data set], https://doi.org/10.16904/envidat.479, 2024. a

Monin, A. S.: The Atmospheric Boundary Layer, Annu. Rev. Fluid Mech., 2, 225–250, https://doi.org/10.1146/annurev.fl.02.010170.001301, 1970. a

Monin, A. S. and Obukhov, A. M.: Basic laws of turbulent mixing in the surface layer of the atmosphere, Transactions of the Geophysical Institute, Academy of Sciences of the USSR, 24, 163–187, 1954. a

Mott, R. and Lehning, M.: Meteorological Modeling of Very High-Resolution Wind Fields and Snow Deposition for Mountains, J. Hydrometeorol., 11, 934–949, https://doi.org/10.1175/2010JHM1216.1, 2010. a

Mott, R., Vionnet, V., and Grünewald, T.: The Seasonal Snow Cover Dynamics: Review on Wind-Driven Coupling Processes, Front. Earth Sci., 6, 197, https://doi.org/10.3389/feart.2018.00197, 2018. a

Nalpanis, P., Hunt, J. C. R., and Barrett, C. F.: Saltating particles over flat beds, J. Fluid Mech., 251, 661–685, https://doi.org/10.1017/S0022112093003568, 1993. a, b, c

Namikas, S. L.: Field measurement and numerical modelling of aeolian mass flux distributions on a sandy beach, Sedimentology, 50, 303–326, https://doi.org/10.1046/j.1365-3091.2003.00556.x, 2003. a, b, c, d, e, f, g, h, i, j, k

Nemoto, M. and Nishimura, K.: Numerical simulation of snow saltation and suspension in a turbulent boundary layer, J. Geophys. Res.-Atmos., 109, D18206, https://doi.org/10.1029/2004JD004657, 2004. a, b, c, d

Niiya, H. and Nishimura, K.: Spatiotemporal Structure of Aeolian Particle Transport on Flat Surface, J. Phys. Soc. Jpn., 86, 054402, https://doi.org/10.7566/JPSJ.86.054402, 2017. a, b

Niiya, H. and Nishimura, K.: Hysteresis and Surface Shear Stresses During Snow-Particle Aeolian Transportation, Bound.-Lay. Meteorol., 183, 447–467, https://doi.org/10.1007/s10546-022-00688-8, 2022. a, b, c, d, e

Nishimura, K. and Hunt, J. C. R.: Saltation and incipient suspension above a flat particle bed below a turbulent boundary layer, J. Fluid Mech., 417, 77–102, https://doi.org/10.1017/S0022112000001014, 2000. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p

Nishimura, K. and Nemoto, M.: Blowing snow at Mizuho station, Antarctica, Philos. T. R. Soc. A, 363, 1647–1662, https://doi.org/10.1098/rsta.2005.1599, 2005. a, b, c, d, e, f, g, h, i, j

Nishimura, K., Yokoyama, C., Ito, Y., Nemoto, M., Naaim-Bouvet, F., Bellot, H., and Fujita, K.: Snow particle speeds in drifting snow: Snow particle speeds in drifting snow, J. Geophys. Res.-Atmos., 119, 9901–9913, https://doi.org/10.1002/2014JD021686, 2014. a, b, c, d

Owen, P. R.: Saltation of uniform grains in air, J. Fluid Mech., 20, 225–242, https://doi.org/10.1017/S0022112064001173, 1964. a, b, c, d, e, f, g, h, i, j, k, l, m

Palm, S. P., Kayetha, V., Yang, Y., and Pauly, R.: Blowing snow sublimation and transport over Antarctica from 11 years of CALIPSO observations, The Cryosphere, 11, 2555–2569, https://doi.org/10.5194/tc-11-2555-2017, 2017. a, b, c, d

Paterna, E., Crivelli, P., and Lehning, M.: Decoupling of mass flux and turbulent wind fluctuations in drifting snow, Geophys. Res. Lett., 43, 4441–4447, https://doi.org/10.1002/2016GL068171, 2016. a, b

Pomeroy, J. W. and Gray, D. M.: Saltation of snow, Water Resour. Res., 26, 1583–1594, https://doi.org/10.1029/WR026i007p01583, 1990. a, b, c, d, e, f, g, h, i, j, k, l

Pomeroy, J. W. and Male, D. H.: Steady-state suspension of snow, J. Hydrol., 136, 275–301, https://doi.org/10.1016/0022-1694(92)90015-N, 1992. a, b

Pomeroy, J. W., Gray, D. M., and Landine, P. G.: The Prairie Blowing Snow Model: characteristics, validation, operation, J. Hydrol., 144, 165–192, https://doi.org/10.1016/0022-1694(93)90171-5, 1993. a

Pomeroy, J. W., Marsh, P., and Gray, D. M.: Application of a distributed blowing snow model to the Arctic, Hydrol. Process., 11, 1451–1464, https://doi.org/10.1002/(SICI)1099-1085(199709)11:11<1451::AID-HYP449>3.0.CO;2-Q, 1997. a

Prandtl, L.: The mechanics of viscous flows, in: Aerodynamic Theory III, edited by Durand, W. F., Springer Berlin, pp. 34–208, https://doi.org/10.1007/978-3-642-91489-8, 1935. a

Raderschall, N., Lehning, M., and Schär, C.: Fine-scale modeling of the boundary layer wind field over steep topography, Water Resour. Res., 44, W09425, https://doi.org/10.1029/2007WR006544, 2008. a

Rasmussen, K. R. and Sørensen, M.: Vertical variation of particle speed and flux density in aeolian saltation: Measurement and modeling, J. Geophys. Res.-Earth, 113, F02S12, https://doi.org/10.1029/2007JF000774, 2008. a

Raupach, M. R.: Saltation layers, vegetation canopies and roughness lengths, in: Aeolian Grain Transport 1, edited by: Barndorff-Nielsen, O. E. and Willetts, B. B., Springer Vienna–New York, pp. 83–96, ISBN 978-3-211-82269-2, https://doi.org/10.1007/978-3-7091-6706-9, 1991. a

Raupach, M. R., Gillette, D. A., and Leys, J. F.: The effect of roughness elements on wind erosion threshold, J. Geophys. Res.-Atmos., 98, 3023–3029, https://doi.org/10.1029/92JD01922, 1993. a

Sato, T., Kosugi, K., and Sato, A.: Saltation-layer structure of drifting snow observed in wind tunnel, Ann. Glaciol., 32, 203–208, https://doi.org/10.3189/172756401781819184, 2001. a, b, c, d

Schiller, L. and Nauman, A. Z.: Über die grundlegenden Berechnungen bei der SchwerSraftaufbereitung, Z. Ver. Dtsch. Ing., 77, 318–320, 1933. a

Schmidt, R.: Estimates of threshold windspeed from particle sizes in blowing snow, Cold Reg. Sci. Technol., 4, 187–193, https://doi.org/10.1016/0165-232X(81)90003-3, 1981. a

Schmidt, R. A.: Vertical profiles of wind speed, snow concentration, and humidity in blowing snow, Bound.-Lay. Meteorol., 23, 223–246, https://doi.org/10.1007/BF00123299, 1982. a

Shao, Y. and Mikami, M.: Heterogeneous Saltation: Theory, Observation and Comparison, Bound.-Lay. Meteorol., 115, 359–379, https://doi.org/10.1007/s10546-004-7089-2, 2005. a

Sharma, V., Parlange, M. B., and Calaf, M.: Perturbations to the Spatial and Temporal Characteristics of the Diurnally-Varying Atmospheric Boundary Layer Due to an Extensive Wind Farm, Bound.-Lay. Meteorol., 162, 255–282, https://doi.org/10.1007/s10546-016-0195-0, 2017. a

Sharma, V., Comola, F., and Lehning, M.: On the suitability of the Thorpe–Mason model for calculating sublimation of saltating snow, The Cryosphere, 12, 3499–3509, https://doi.org/10.5194/tc-12-3499-2018, 2018. a, b, c, d

Sharma, V., Gerber, F., and Lehning, M.: Introducing CRYOWRF v1.0: multiscale atmospheric flow simulations with advanced snow cover modelling, Geosci. Model Dev., 16, 719–749, https://doi.org/10.5194/gmd-16-719-2023, 2023. a, b, c, d, e, f, g, h, i, j, k, l

Sigmund, A., Dujardin, J., Comola, F., Sharma, V., Huwald, H., Melo, D. B., Hirasawa, N., Nishimura, K., and Lehning, M.: Evidence of Strong Flux Underestimation by Bulk Parametrizations During Drifting and Blowing Snow, Bound.-Lay. Meteorol., 182, 119–146, https://doi.org/10.1007/s10546-021-00653-x, 2022. a, b

Skamarock, W. C., Klemp, B., Dudhia, J., Gill, O., Liu, Z., Berner, J., Wang, W., Powers, G., Duda, G., Barker, D., and Huang, X.-y.: A Description of the Advanced Research WRF Model Version 4.1, Tech. Rep. No. NCAR/TN-556+STR, National Center for Atmospheric Research, Boulder, Colorado, https://doi.org/10.5065/1dfh-6p97, 2019. a, b

Sommer, C. G., Wever, N., Fierz, C., and Lehning, M.: Investigation of a wind-packing event in Queen Maud Land, Antarctica, The Cryosphere, 12, 2923–2939, https://doi.org/10.5194/tc-12-2923-2018, 2018. a

Sørensen, M.: An analytic model of wind-blown sand transport, in: Aeolian Grain Transport 1, edited by: Barndorff-Nielsen, O. E. and Willetts, B. B., Springer Vienna–New York, ISBN 978-3-211-82269-2, https://doi.org/10.1007/978-3-7091-6706-9, pp. 67–81, 1991. a, b, c, d, e

Sørensen, M.: On the rate of aeolian sand transport, Geomorphology, 59, 53–62, https://doi.org/10.1016/j.geomorph.2003.09.005, 2004. a, b, c, d, e, f

Stull, R. B.: An Introduction to Boundary Layer Meteorology, Kluwer Academic Publishers, Dordrecht, The Netherlands, ISBN-13 978-90-277-2769-5, https://doi.org/10.1007/978-94-009-3027-8, 1988. a

Sugiura, K., Nishimura, K., Maeno, N., and Kimura, T.: Measurements of snow mass flux and transport rate at different particle diameters in drifting snow, Cold Reg. Sci. Technol., 27, 83–89, https://doi.org/10.1016/S0165-232X(98)00002-0, 1998. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w

Takeuchi, M.: Vertical profile and Horizontal Increase of Drift-Snow Transport, J. Glaciol., 26, 481–492, https://doi.org/10.3189/S0022143000010996, 1980. a

Thorpe, A. D. and Mason, B. J.: The evaporation of ice spheres and ice crystals, Brit. J. Appl. Phys., 17, 541–548, https://doi.org/10.1088/0508-3443/17/4/316, 1966. a

Ungar, J. E. and Haff, P. K.: Steady state saltation in air, Sedimentology, 34, 289–299, https://doi.org/10.1111/j.1365-3091.1987.tb00778.x, 1987. a, b, c, d, e, f, g

van den Broeke, M., König-Langlo, G., Picard, G., Kuipers Munneke, P., and Lenaerts, J.: Surface energy balance, melt and sublimation at Neumayer Station, East Antarctica, Antarct. Sci., 22, 87–96, https://doi.org/10.1017/S0954102009990538, 2010. a

van Wessem, J. M., van de Berg, W. J., Noël, B. P. Y., van Meijgaard, E., Amory, C., Birnbaum, G., Jakobs, C. L., Krüger, K., Lenaerts, J. T. M., Lhermitte, S., Ligtenberg, S. R. M., Medley, B., Reijmer, C. H., van Tricht, K., Trusel, L. D., van Ulft, L. H., Wouters, B., Wuite, J., and van den Broeke, M. R.: Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 2: Antarctica (1979–2016), The Cryosphere, 12, 1479–1498, https://doi.org/10.5194/tc-12-1479-2018, 2018. a, b, c, d, e, f, g, h

Vionnet, V., Martin, E., Masson, V., Guyomarc'h, G., Naaim-Bouvet, F., Prokop, A., Durand, Y., and Lac, C.: Simulation of wind-induced snow transport and sublimation in alpine terrain using a fully coupled snowpack/atmosphere model, The Cryosphere, 8, 395–415, https://doi.org/10.5194/tc-8-395-2014, 2014. a, b, c, d, e, f, g, h, i, j, k, l

Walter, B., Horender, S., Voegeli, C., and Lehning, M.: Experimental assessment of Owen's second hypothesis on surface shear stress induced by a fluid during sediment saltation, Geophys. Res. Lett., 41, 6298–6305, https://doi.org/10.1002/2014GL061069, 2014. a, b

Werner, B. T.: A Steady-State Model of Wind-Blown Sand Transport, J. Geol., 98, 1–17, 1990. a

White, P. W.: Part IV: Physical processes, in: IFS Documentation CY23R4, vol. 4, European Centre for Medium-Range Weather Forecasts (ECMWF), https://doi.org/10.21957/02054f0fbf, 2003. a

Yu, H., Li, G., Walter, B., Lehning, M., Zhang, J., and Huang, N.: Wind conditions for snow cornice formation in a wind tunnel, The Cryosphere, 17, 639–651, https://doi.org/10.5194/tc-17-639-2023, 2023. a