Articles | Volume 12, issue 12
Research article
10 Dec 2018
Research article |  | 10 Dec 2018

A simulation of a large-scale drifting snowstorm in the turbulent boundary layer

Zhengshi Wang and Shuming Jia

Drifting snowstorms are an important aeolian process that reshape alpine glaciers and polar ice shelves, and they may also affect the climate system and hydrological cycle since flying snow particles exchange considerable mass and energy with air flow. Prior studies have rarely considered full-scale drifting snowstorms in the turbulent boundary layer; thus, the transportation feature of snow flow higher in the air and its contribution are largely unknown. In this study, a large-eddy simulation is combined with a subgrid-scale velocity model to simulate the atmospheric turbulent boundary layer, and a Lagrangian particle tracking method is adopted to track the trajectories of snow particles. A drifting snowstorm that is hundreds of meters in depth and exhibits obvious spatial structures is produced. The snow transport flux profile at high altitude, previously not observed, is quite different from that near the surface; thus, the extrapolated transport flux profile may largely underestimate the total transport flux. At the same time, the development of a drifting snowstorm involves three typical stages, rapid growth, gentle growth, and equilibrium, in which large-scale updrafts and subgrid-scale fluctuating velocities basically dominate the first and second stages, respectively. This research provides an effective way to gain an insight into natural drifting snowstorms.

1 Introduction

Snow, one type of solid precipitation, is an important source of material to mountain glaciers and polar ice sheets, which are widespread throughout high and cold regions (Chang et al., 2016; Gordon and Taylor, 2009; Lehning et al., 2008). A common natural phenomenon over snow cover is the drifting snowstorm, which occurs when the wind speed exceeds a critical value (Doorschot et al., 2004; Li and Pomeroy, 1997; Sturm and Stuefer, 2013). Drifting snow can entrain loose snow particles on the bed into the air, which may be further transported to high altitude by turbulent eddies (King, 1990; Mann et al., 2000; Nemoto and Nishimura, 2004). Drifting snow clouds can typically range in thickness from tens to thousands of meters (Mahesh et al., 2003; Palm et al., 2011), which may not only affect people's daily life by reducing the visibility and producing local accumulation (Gordon and Taylor, 2009; Mohamed et al., 1998), but can also influence the global climate system evolution by changing the mass and energy balance of ice shelves (Cess and Yagai, 1991; Hanesiak and Wang, 2005; Hinzman et al., 2005; Lenaerts and Broeke, 2012).

Several field experiments on drifting snowstorms have been performed (Bintanja, 2001; Budd, 1966; Dingle and Radok, 1961; Doorschot et al., 2004; Gallée et al., 2013; Gordon and Taylor, 2009; Guyomarch et al., 2014; Kobayashi, 1978; Mann et al., 2000; Nishimura and Nemoto, 2005; Nishimura et al., 2015; Pomeroy and Gray, 1990; Sbuhei, 1985; Schmidt, 1982; Sturm and Stuefer, 2013) since the middle of the last century. However, the measurements are commonly conducted near the surface; thus, drifting snow features at high altitude are unknown, and the impacts of these features are difficult to assess. A thorough investigation documenting the evolution process and structure of a full-scale drifting snowstorm is essential to understand this natural phenomenon and assess its impacts.

Drifting snow models, however, offer a panoramic view of the evolution process of drifting snow and thus have become one of the most useful research approaches. Many continuum medium models of drifting snow (Bintanja, 2000; Déry and Yau, 1999; Schneiderbauer and Prokop, 2011; Uematsu et al., 1991; Vionnet et al., 2014) have advanced the knowledge of natural drifting snow to a great extent. However, a particle-tracking drifting snow model is still needed since the particle characteristics and its motion require further investigation. Although a series of particle tracking models (Huang et al., 2016; Huang and Shi, 2017; Huang and Wang, 2015, 2016; Nemoto and Nishimura, 2004; Zhang and Huang, 2008; Zwaaftink et al., 2014) have been established, these models have generally focused on the grain–bed interactions and particle motions near the surface. Thus, a drifting snow model aimed at producing a large-scale drifting snowstorm in a turbulent boundary layer deserves further exploration.

In this study, a drifting snow model in the atmospheric boundary layer that focuses on a full-scale drifting snowstorm is established. The wind field is solved using a large-eddy simulation for the purpose of generating a turbulent atmospheric boundary layer. A subgrid-scale (SGS) velocity is also considered to include the diffusive effect of small-scale turbulence. Finally, particle motion is calculated using a Lagrangian particle tracking method. The large-scale drifting snowstorm is produced under the actions of large-scale turbulent structures combined with a steady-state snow saltation boundary condition for particles, and its spatial structures and transport features are analyzed.

2 Model and methods

2.1 Simulation of a turbulent atmospheric boundary layer

The mesoscale atmosphere prediction pattern ARPS (Advanced Regional Prediction System, version 5.3.3) is adopted to simulate the turbulent atmospheric boundary layer, in which the filtered three-dimensional compressible non-hydrostatic Navier–Stokes equation is solved (Xue et al., 2001):


where “” represents variables that are filtered and the filtering scale is Δ̃=Δx1Δx2Δx31/3, in which Δxi is the grid spacing along the streamwise (i=1), spanwise (i=2), and vertical directions (i=3). ρ=p1-qv/ε+qv1+qv/RdT is the air density, in which p, qv, Rd, and T are the pressure, specific humidity, gas constant (287.0 J kg−1 K−1), and temperature of the air, respectively, and ε=0.622 is a constant. ui is the instantaneous wind speed component, and xi is the position coordinate. t is time, δij is the Kronecker delta, B=-gρ/ρ is the buoyancy caused by the air density perturbation ρ, and g is the acceleration due to gravity. p*=p-α(ρu) contains the pressure perturbation term and damping term, where α=0.5 is the damping coefficient and  is the divergence. The subgrid stress τij can be expressed as (Smagorinsky, 1963)

(3) τ i j = - 2 ν t S ̃ i j = - 2 C s Δ ̃ 2 | S ̃ | S ̃ i j ,

where S̃ij=0.5ũi/xj+ũj/xi is the strain rate tensor and |S̃|=2S̃ijS̃ij, Cs is the Smagorinsky coefficient that is determined locally by the dynamic Lagrangian model (Meneveau et al., 1996).

2.2 Governing equation of particle motion

The trajectory of each snow particle is calculated using a Lagrangian particle tracking method. Since a snow particle has is almost 103 times more dense than air, airborne particles are assumed to process only gravity and fluid drag forces, and the governing equations of particle motion can be expressed as (Dupont et al., 2013; Huang and Wang, 2016; Vinkovic et al., 2006)


where xpi and upi are the position coordinate and velocity of the snow particle, respectively. mp is the mass of the solid particle, Vr is the relative speed between the snow particle and air, and Tp=ρpdp2/18ρν is the particle relaxation time, where ρp is the particle density (900 kg m−3), dp is the particle diameter and ν=1.5×10-5 m2 s−1 is the kinematic viscosity of air. f(Rep) can be expressed as (Clift et al., 1978)

(6) f R e p = 1 R e p < 1 1 + 0.15 R e p 0.687 R e p 1 ,

where Rep=Vrd/ν is the particle Reynolds number.

Considering the large grid spacing in simulating an atmospheric boundary layer (for which the information about turbulent vortices smaller than the grid size is missing), the SGS velocity is also included and attached to the particle. Namely, the local relative is expressed as Vri=ũixp-upi+ui, in which ũixp is the resolved large-scale wind speed at the particle's position and is determined by the resolved wind speeds of surrounding grid points through the linear interpolation algorithm. The SGS velocity can be calculated from the SGS stochastic model of Vinkovic et al. (2006):

(7) d u i = - 1 T L + 1 2 k ̃ d k ̃ d t u i d t + 4 k ̃ 3 T L d η i ( t ) ,

where TL=4k̃/3C0ε̃ is the Lagrangian correlation timescale. Here, C0=2.1 is the Lagrangian constant, ε̃=Cεk̃3/2/Δ̃ is the subgrid turbulence dissipation rate, Cε=0.41 is a constant, and dηi is the increment of a vector-valued Wiener process with zero mean and variance dt. k̃ is the subgrid turbulent kinetic energy and can be obtained from the transport equation (Deardorff, 1980)

(8) k ̃ t + u ̃ j k ̃ x j = ν t 3 g θ 0 θ ̃ x 3 + 2 ν t S ̃ i j 2 + 2 x j ν t k ̃ x j + ε ̃ ,

where θ is the potential temperature and θ0 is the surface potential temperature.

2.3 Initial conditions of snow particles

To generate a large-scale drifting snowstorm, a steady-state snow saltation condition is set as the bottom boundary condition for particles. During drifting snow events, the sum of residual fluid shear stress τf and particle-borne shear stress τp should be equal to the total shear stress τ; thus, the particle-borne stress can be expressed as

(9) τ p = τ - τ f .

Here, the residual fluid shear stress τf is set to be the threshold shear stress τtf of drifting snow, which can be read as (Clifton et al., 2006)

(10) τ tf = A 2 g d p ρ p - ρ ,

in which A=0.2 is a constant, and dp is the mean diameter of the snow particles.

At the same time, the particle-borne shear stress at the surface can be calculated from the particle trajectories as (Nemoto and Nishimura, 2004)

(11) τ p = i = 1 n m i u p i - i = 1 n m i u p i ,

where mi is the mass of particle and upi and upi are the horizontal speeds of impact and lift-off particles, respectively. n and n are the particle numbers per unit area in unit time of impact and lift-off grains, respectively, which should be equivalent in steady-state saltation. Thus, the number of lift-off particles per unit area is

(12) n = n = τ p m i 1 - e h u p i ,

in which 〈 〉 indicates the overall average, and eh is the horizontal restitution coefficient of snow particles. According to Sugiura and Maeno (2000), the mean horizontal restitution coefficient can be expressed as

(13) e h = 0.48 θ i 0.01 v i 1.27 m s - 1 0.48 v i 1.27 - log v i 1.27 θ i 0.01 v i > 1.27 m s - 1 ,

where θi and vi are the impact velocity and angle, respectively. Here, θi has a mean value of approximately 10 (Sugiura and Maeno, 2000), and vi is set to be the threshold of impact velocity. Considering the steady-state saltation condition (one impact particle generates one ejecta on average), vi is determined by setting the ejection number ne=0.51vi0.6θi0.16 equal to 1. In this way, the mean horizontal velocity of impact particles can be obtained through upi=vicosθi.

Then, the velocities of lift-off particles can be obtained from the restitution coefficient of snow. The horizontal restitution coefficient obeys the normal distribution with a mean value given in Eq. (13) and a standard variance as follows (Sugiura and Maeno, 2000):

(14) σ 2 = 0.07 θ i - 0.06 v i 0.52 m s - 1 0.07 v i 0.52 - log v i 0.52 θ i - 0.06 v i > 0.52 m s - 1 .

Conversely, the vertical restitution coefficient can be described by a two parameter gamma function (see Eq. 17), in which the parameters α and β can be expressed as (Sugiura and Maeno, 2000)



In this condition, if some of the snow particles within the saltation layer are transported to higher in the air by turbulent vortexes (the saltation layer becomes undersaturated), more particles will lift off from the surface to replenish the saltation layer until a saturated state is reached.

2.4 Simulation details

The computational domain is 1000 m × 500 m × 1000 m, with a uniform horizontal grid size of 5 m adopted to solve finer vortex structure in the atmospheric boundary layer. The mean grid size in the vertical direction is 20 m, with a grid refinement algorithm adopted near the surface (the finest grid size is 1 m). Periodic boundaries are used along streamwise and spanwise dimensions, and the bottom is set as a grid wall. The top is set as an open radiation boundary with a Rayleigh damping layer that is 250 m in depth.

The atmosphere is neutral with an initial potential temperature of 300 K and an initial relative humidity of 90 %. The initial wind profile is logarithmic with a surface roughness of 0.1 m (Doorschot et al., 2004). Atmospheric turbulence is induced by random initial potential temperature perturbations at the first-level grid level with a maximum magnitude of 0.5 K and is sustained by a constant heat flux at the bottom. The constant heat flux is 50 W m−2 according to the observation of Pomeroy and Essery (1999). The evolution time for a turbulent boundary layer is 5 times the large-eddy turnover time t*(H/u*, where H is the boundary layer depth and u* is the friction velocity). Actually, this condition corresponds to an “intermediate” turbulent boundary layer that is dominated by wind shear force (Moeng and Sullivan, 1994). Thus, the structures of the drifting snowstorm should not be affected by the changing surface heat flux significantly if the surface heat flux is small. Further simulations with different values of surface heat flux (<100 W m−2) also prove this point.

For particles, periodic boundary conditions are also used at lateral boundaries, and a rebound boundary condition without energy loss is adopted at the model top. The bottom boundary condition for particles is given in Sect. 2.3 and is updated every 0.5 s. Additionally, each particle represents one particle parcel for the purpose of reducing computational complexity. In this simulation, each particle parcel contains 107 snow particles. The large time step and small time step (acoustic wave integral) for the wind field calculation are 0.1 and 0.02 s, respectively, and the particle time step is determined by the minimum of particle relaxation time.

The size distribution of lift-off particles in drifting snow can be described well by the two-parameter gamma function (Budd, 1966; Gordon and Taylor, 2009; Nishimura and Nemoto, 2005; Schmidt, 1982)

(17) f d p = d p α p - 1 β p α p Γ α p exp - β p d p ,

where αp and βp are the shape and scale parameters of the distribution, respectively. In this simulation, the diameters of lift-off snow particles are given randomly from a gamma function with the parameters of αp=4 and βp=50, as shown in Fig. 1, which is also consistent with observed particle size distributions (Nishimura and Nemoto, 2005; Schmidt, 1982).

Figure 1Size distribution of lift-off snow particles in this simulation.


3 Results and discussions

3.1 Model validation

When drifting snow occurs in the atmospheric boundary layer, updrafts and turbulence fluctuations can send snow particles to high altitude, forming a fully developed drifting snowstorm. Figure 2 shows the drifting snowstorm in the atmospheric boundary layer at different moments, in which the friction velocity is u*=0.29 m s−1 and dark spots represent snow particles. It can be seen that drifting snowstorms experience an evolution process from near the surface to high altitudes, in which particle concentration decreases with increasing height. The high concentrations of drifting snow cloud are generally below 500 m, though snow particles may reach up to approximately 800 m under this condition. This is also consistent with observations (Mahesh et al., 2003; Palm et al., 2011).

Figure 2Drifting snowstorm at different moments under the friction velocity of 0.29 m s−1.


Since a drifting snowstorm exhibits a different structure from bottom to top, the evolution of particle number density profile in the drifting snowstorm is shown in Fig. 3, which is also compared with measurements of Mann et al. (2000). From this figure, the thickness of the drifting snow layer obviously increases with time and almost approaches its steady state after 1200 s. At the same time, the particle number density basically decreases with height, which is consistent with the measurements of Mann et al. (2000) at various friction velocities. The predicted particle number density at the surface is much larger than at higher altitude and observations, mainly because the saltating particles are also included.

Generally, smaller particles are more likely to be transported higher in the air. Figure 4 shows the variation in modeled average particle diameter versus height, which is also compared with various field measurements (Nishimura and Nemoto, 2005; Schmidt, 1982). Similar to the field observations, the average particle size basically decreases with height at lower altitude but is almost constant above 1 m. The average particle diameter is approximately 75 µm, ranging from 1 m to hundreds of meters in height, which is also consistent with the measurements of Nishimura and Nemoto (2005).

Figure 3Evolution of particle number density under the friction velocities (a) 0.29 m s−1 and (b) 0.51 m s−1.


Then, the particle size distributions at various heights are also compared with experiment results. As shown in Fig. 5, the heights are 0.05, 0.5, and 1 m. The modeled particle size distributions at various heights are consistent with the measurements (Nishimura and Nemoto, 2005; Schmidt, 1982). Therefore, the established model is able to produce a large-scale drifting snowstorm.

In addition, it can be seen that the proportion of particles below 100 µm in diameter at 0.05 m is smaller than that of the experimental result. The reason could be that midair collisions, occurring frequently within the high-concentration saltating snow cloud at the near surface, play an important role in conveying larger particles to higher altitude (Carneiro et al., 2013). However, the midair collision mechanism is beyond the scope of the current study.

3.2 Snow transport flux

The snow transport flux is of great importance to predict the mass and energy balances of ice sheets. The total transport flux can be obtained from vertical integration of the snow transport flux profile.

Figure 4Variation in average particle diameter versus height.


The profiles of snow transport rate, per unit area, per unit time, under various friction velocities are shown in Fig. 6a. It can be seen that the transport flux undergoes a sharp decrease with height at a lower altitude (e.g., below 1.0 m); however, the transport flux tends to decrease rather gently until almost the top of the drifting snowstorm, as shown in Fig. 6b, probably due to the large-scale turbulent motion and increasing wind speed with height. In other words, the suspension flux of drifting snow at higher altitudes, previously not observed, may be much larger than we previously thought. The mean horizontal wind speed profiles of the fully developed turbulent boundary layer under various friction velocities are shown in Fig. 7. The horizontal wind speed increases with height and changes into a constant above the boundary layer. The rapid decrease in the snow transport flux occurs at about the top of the turbulent boundary layer, mainly because turbulence becomes weaker above this height and fewer particles can be transported to a higher altitude.

In addition, the transition of snow transport flux profile at about 1 m should be mainly caused by the different motion states of particles with different particle sizes, as shown in Fig. 4. Above the critical height, particles generally follow the turbulent flow in the state of suspension because their gravities and relaxation times are small enough. However, plenty of larger particles at the near surface make the particles' velocity differ from the wind speed since particle inertia plays an important role.

In previous studies, only the transport fluxes at the near surface are commonly measured (Mann et al., 2000; Nishimura and Nemoto, 2005; Schmidt, 1982, 1984; Tabler, 1991); thus, the features of the entire transport flux profile are largely unclear, which may result in considerable uncertainties about the total transport flux. The proportions of suspension flux above a given height hc (referred as Qc) to the total suspension flux Qs are shown in Fig. 7, in which snow particles below 0.1 m are not calculated (Mann et al., 2000).

From Fig. 8a, the contribution of Qc to the total suspension flux is non-negligible under various hc values, the proportion of Qc when hc=100 m to the total suspension flux has exceeded 30 % when the friction velocity is 0.46 m s−1. At the same time, the proportion of Qc to the total suspension flux increases with friction velocity but decreases with increasing hc. From Fig. 8b, it can be seen that the proportion of Qc to the total suspension flux is only slightly affected by the surface heat flux, which indicates that the structures of drifting snowstorm are not sensitive to the surface heat flux under this condition. The influence of surface heat flux is also weakened by the increasing friction velocity, mainly because larger friction velocity results in stronger turbulence under the actions of wind shear.

In this way, not only the snow transport flux, but also the sublimation of suspended snow particles should be reevaluated because the sublimation rate of snow particles higher in the air may be much larger than near the surface due to the lower air humidity and greater wind speed at higher altitude (Mann et al., 2000; Nishimura and Nemoto, 2005; Schmidt, 1982, 1984; Tabler, 1991).

3.3 Structures in a drifting snowstorm

In a drifting snowstorm, particles aggregate locally and produce special spatial structures (as shown in Fig. 2). These structures should be directly related to the turbulence structures present in the atmospheric boundary layer. Drifting snowstorms without atmospheric turbulence are shown in Fig. 9. This simulation is achieved by replacing the resolved wind speed at a particle's position ũixp with a given value obtained from the standard logarithmic profile, and the other model settings and simulation procedures stay the same with other simulations. In this way, the effect of large-scale turbulent structures on the development of the drifting snowstorm vanishes. Compared with Fig. 2, drifting snow particles mainly travel at the near surface with a uniform spatial distribution when atmospheric turbulence is not included.

Figure 5Particle size distribution at various heights.


It is known that snow particles will become suspended if the local vertical wind speed exceeds the terminal velocity of a particle. In a turbulent atmospheric boundary layer, there is a large number of turbulent structures with different scales and shapes. The vertical wind speed component of large-scale turbulence (namely, updraft) plays an important role in carrying snow particles to high altitude, while small-scale turbulence (e.g., the SGS fluctuating velocity) tends to spread particles from high concentration zones to low concentration zones. As shown in Fig. 10a, at the initial period of a drifting snowstorm, the structures in the drifting snowstorm are consistent with large-scale updrafts, and snow particles are mainly located in the updraft. With the further development of the drifting snowstorm, as shown in Fig. 10b, more snow particles are scattered around the updraft bubbles, although high-concentration particle clouds are still in the wind bubbles. When a drifting snowstorm approaches its saturated state, snow particle clouds are almost connected together with numerous high-concentration zones inside.

Figure 6Variations in snow transport flux versus height.


Figure 7Horizontal wind speed profiles of the fully developed turbulent boundary layer under various friction velocities.


The evolution of the depth of a drifting snowstorm can be divided into three typical stages. In sequence, these phases are the rapid growth phase, the gentle growth stage, and an equilibrium state, as shown in Fig. 11. Here, the depth of a drifting snowstorm refers to the average height of the topmost particle during this period (100 s). The rapid growth stage is mainly driven by large-scale turbulent motion, while the turbulent diffusion by the SGS fluctuating velocity is the main contributor to the gentle growth stage. The duration of the second stage decreases with increasing friction velocity, which mainly comes from the stronger turbulent diffusion under larger friction velocities.

Figure 8Proportion of suspension flux above hc to the total suspension flux under (a) various friction velocities and (b) various surface heat fluxes Qs.


Figure 9Drifting snowstorm without atmospheric turbulence under a friction velocity of 0.35 m s−1.


Figure 10Evolution of the drifting snowstorm and vertical wind speed bubbles under a friction velocity of 0.35 m s−1, and wind bubbles are the iso-surface of vertical wind speed with a value of 1.0 m s−1 (corresponding to the critical wind speed at which the particle of mean particle size becomes a suspended particle since the maximum diameter of suspended particles is found to be approximately equal to the mean particle size of the lift-off particles).


At the same time, the time required for the drifting snowstorm to reach its maximum thickness decreases with friction velocity, ranging from about 1200 s to approximately 600 s when the friction velocity increases from 0.29 to 0.46 m s−1. The thicknesses of saturated drifting snowstorms are almost constant with a value of approximately 900 m under different friction velocities, probably because the boundary layer depth as well as the surface heat flux are unchanged. Higher domain heights are also tested with the same model settings, and the thickness of the drifting snow seems basically unchanged. Drifting snowstorms with a difference in thickness may be achieved by changing the initial state of the air and surface heat flux. Thus, the final thickness of a drifting snowstorm should be largely dependent on the maximum height of atmospheric turbulence.

Figure 11Time evolutions of the thickness of the drifting snowstorm under various friction velocities.


4 Conclusions

In this work, large-scale drifting snowstorms are simulated in a large-eddy simulation combined with a particle tracking model that includes subgrid-scale velocity fluctuations. A typical drifting snowstorm of several hundred meters in depth is generated, and the structure of the particle cloud with different concentrations is also produced. The transport flux profile has obviously different slopes near the surface compared to higher altitudes; that is, transport flux at the near surface decreases with height sharply, but decreases more gently at higher altitude. Previous studies may largely underestimate the total transport during drifting snowstorms.

At the same time, the evolution of the thickness of a drifting snowstorm generally contains three stages. Drifting snowstorm development generally begins with a rapid growth stage driven by large-scale atmospheric turbulent motions, followed by a gentle growth stage driven by the SGS fluctuating wind speed, before reaching an equilibrium stage when the drifting snow approaches a saturated state. The second stage becomes shorter with increasing friction velocity, mainly because stronger turbulence under higher friction velocity enhances the turbulent diffusion of particles.

Data availability

The data used in this publication are available upon request from the corresponding author. The source code of ARPS is publically available online at the Center for Analysis and Prediction of Storms at the University of Oklahoma (, Xue et al., 2001).

Author contributions

ZW and SJ designed the research. ZW implemented the computational model and analysis tools, carried out the simulations, and collected the data. ZW and SJ wrote the paper. All authors reviewed the paper.

Competing interests

The authors declare that they have no conflict of interest.


This work is supported by the CARDC Fundamental and Frontier Technology Research Fund (FFTRF-2017-08, FFTRF-2017-09), the National Natural Science Foundation of China (11772143, 41371034), and the National Key Research and Development Program of China (2016YFC0500900).

Edited by: Valentina Radic
Reviewed by: two anonymous referees


Bintanja, R.: Snowdrift suspension and atmospheric turbulence. Part I: Theoretical background and model description, Bound.-Lay. Meteorol., 95, 343–368, 2000. 

Bintanja, R.: Characteristics of snowdrift over a bare ice surface in Antarctica, J. Geophys. Res.-Atmos., 106, 9653–9659, 2001. 

Budd, W. F.: The Byrd snow drift project : outline and basic results, American Geophysical Union, Washington, D.C., 71–134, 1966. 

Carneiro, M. V., Araújo, N. A., Pähtz, T., and Herrmann, H. J.: Midair collisions enhance saltation, Phys. Rev. Lett., 111, 058001,, 2013. 

Cess, R. D. and Yagai, I.: Interpretation of Snow-Climate Feedback as Produced by 17 General Circulation Models, Science, 253, 888–892, 1991. 

Chang, A. T. C., Foster, J. L., and Hall, D. K.: Nimbus-7 SMMR Derived Global Snow Cover Parameters, Ann. Glaciol., 9, 39–44, 2016. 

Clift, R., Grace, J. R., and Weber, M. E.: Bubbles, drops, and particles, Academic Press, New York, 1978. 

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

Deardorff, J. W.: Stratocumulus-capped mixed layers derived from a three-dimensional model, Bound.-Lay. Meteorol., 18, 495–527, 1980. 

Déry, S. J. and Yau, M. K.: A Bulk Blowing Snow Model, Bound.-Lay. Meteorol., 93, 237–251, 1999. 

Dingle, W. R. J. and Radok, U.: Antarctic snow drift and mass transport, Int. Assoc. Sci. Hydrol. Publ., 55, 77–81, 1961. 

Doorschot, J. J. J., Lehning, M., and Vrouwe, A.: Field measurements of snow-drift threshold and mass fluxes, and related model simulations, Bound.-Lay. Meteorol., 113, 347–368, 2004. 

Dupont, S., Bergametti, G., Marticorena, B., and Simoëns, S.: Modeling saltation intermittency, J. Geophys. Res.-Atmos., 118, 7109–7128, 2013. 

Gallée, H., Trouvilliez, A., Agosta, C., Genthon, C., Favier, V., and Naaim-Bouvet, F.: Transport of Snow by the Wind: A Comparison Between Observations in Adélie Land, Antarctica, and Simulations Made with the Regional Climate Model MAR, Bound.-Lay. Meteorol., 146, 133–147, 2013. 

Gordon, M. and Taylor, P. A.: Measurements of blowing snow, Part I: Particle shape, size distribution, velocity, and number flux at Churchill, Manitoba, Canada, Cold Reg. Sci. Technol., 55, 63–74, 2009. 

Guyomarch, G., Goetz, D., Vionnet, V., Naaimbouvet, F., and Deschatres, M.: Observation of Blowing Snow Events and Associated Avalanche Occurrences, in: Preceedings, International Snow Science Workshop, Banff, 2014. 

Hanesiak, J. M. and Wang, X. L.: Adverse-Weather Trends in the Canadian Arctic, J. Climate, 18, 3140–3156, 2005. 

Hinzman, L. D., Bettez, N. D., Bolton, W. R., Chapin, F. S., Dyurgerov, M. B., Fastie, C. L., Griffith, B., Hollister, R. D., Hope, A., and Huntington, H. P.: Evidence and Implications of Recent Climate Change in Northern Alaska and Other Arctic Regions, Climatic Change, 72, 251–298, 2005. 

Huang, N. and Shi, G.: The significance of vertical moisture diffusion on drifting snow sublimation near snow surface, The Cryosphere, 11, 3011–3021,, 2017. 

Huang, N. and Wang, Z.: A 3-D simulation of drifting snow in the turbulent boundary layer, The Cryosphere Discuss., 9, 301–331,, 2015. 

Huang, N. and Wang, Z. S.: The formation of snow streamers in the turbulent atmosphere boundary layer, Aeolian Res., 23, 1–10, 2016. 

Huang, N., Dai, X., and Zhang, J.: The impacts of moisture transport on drifting snow sublimation in the saltation layer, Atmos. Chem. Phys., 16, 7523–7529,, 2016. 

King, J. C.: Some measurements of turbulence over an antarctic ice shelf, Q. J. Roy. Meteorol. Soc., 116, 379–400, 1990. 

Kobayashi, S.: Snow Transport by Katabatic Winds in Mizuho Camp Area, East Antarctica, J. Meteorol. Soc. Jpn., 56, 130–139, 1978. 

Lehning, M., Löwe, H., Ryser, M., and Raderschall, N.: Inhomogeneous precipitation distribution and snow transport in steep terrain, Water Resour. Res., 44, 278–284, 2008. 

Lenaerts, J. T. M. and Broeke, M. R. V. D.: Modeling drifting snow in Antarctica with a regional climate model: 2. Results, J. Geophys. Res.-Atmos., 117, D05109,, 2012. 

Li, L. and Pomeroy, J. W.: Estimates of Threshold Wind Speeds for Snow Transport Using Meteorological Data, J. Appl. Meteorol., 36, 205–213, 1997. 

Mahesh, A., Eager, R., Campbell, J. R., and Spinhirne, J. D.: Observations of blowing snow at the South Pole, J. Geophys. Res.-Atmos., 108, 4707,, 2003. 

Mann, G. W., Anderson, P. S., and Mobbs, S. D.: Profile measurements of blowing snow at Halley, Antarctica, J. Geophys. Res.-Atmos., 105, 24491–24508, 2000. 

Meneveau, C., Lund, T. S., and Cabot, W. H.: A Lagrangian dynamic subgrid-scale model of turbulence, J. Fluid Mech., 319, 353–385, 1996. 

Moeng, C. H. and Sullivan, P. P.: A Comparison of Shear- and Buoyancy-Driven Planetary Boundary Layer Flows, J. Atmos. Sci., 51, 999–1022, 1994. 

Mohamed, N., Florence, N. B., and Hugo, M.: Numerical simulation of drifting snow: erosion and deposition models, Ann. Glaciol., 26, 191–196, 1998. 

Nemoto, M. and Nishimura, K.: Numerical simulation of snow saltation and suspension in a turbulent boundary layer, J. Geophys. Res.-Atmos., 109, D18206,, 2004. 

Nishimura, K. and Nemoto, M.: Blowing snow at Mizuho station, Antarctica, Philos. Trans., 363, 1647–1662, 2005. 

Nishimura, K., Yokoyama, C., Ito, Y., Nemoto, M., Naaim-Bouvet, F., Bellot, H., and Fujita, K.: Snow particle speeds in drifting snow, J. Geophys. Res.-Atmos., 119, 9901–9913, 2015. 

Palm, S. P., Yang, Y., Spinhirne, J. D., and Marshak, A.: Satellite remote sensing of blowing snow properties over Antarctica, J. Geophys. Res.-Atmos., 116, D16123,, 2011. 

Pomeroy, J. W. and Essery, R. L. H.: Turbulent fluxes during blowing snow: field tests of model sublimation predictions, Hydrol. Process., 13, 2963–2975, 1999. 

Pomeroy, J. W. and Gray, D. M.: Saltation of snow, Water Resour. Res., 26, 1583–1594, 1990. 

Sbuhei, T.: Characteristics of Drifting Snow at Mizuho Station, Antarctica, Ann. Glaciol., 6, 71–75, 1985. 

Schmidt, R. A.: Vertical profiles of wind speed, snow concentration, and humidity in blowing snow, Bound.-Lay. Meteorol., 23, 223–246, 1982. 

Schmidt, R. A.: Transport rate of drifting snow and the mean wind speed profile, Bound.-Lay. Meteorol., 34, 213–241, 1984. 

Schneiderbauer, S. and Prokop, A.: The atmospheric snow-transport model: SnowDrift3D, J. Glaciol., 57, 526–542, 2011. 

Smagorinsky, J.: General Circulation Experiments with the Primitive Equations, Mon. Weather Rev., 91, 99–164, 1963. 

Sturm, M. and Stuefer, S.: Wind-blown flux rates derived from drifts at arctic snow fences, J. Glaciol., 59, 21–34, 2013. 

Sugiura, K. and Maeno, N.: Wind-Tunnel Measurements Of Restitution Coefficients And Ejection Number Of Snow Particles In Drifting Snow: Determination Of Splash Functions, Bound.-Lay. Meteorol., 95, 123–143, 2000. 

Tabler, R. D.: Snow transport as a function of wind speed and height, in: Proceedings, Cold Regions Sixth International Specialty Conference TCCP/ASCE, Cold Regions Engineering, 26–28 February 1991, West Lebanon, NH, 729–738, 1991. 

Uematsu, T., Nakata, T., Takeuchi, K., Arisawa, Y., and Kaneda, Y.: Three-dimensional numerical simulation of snowdrift, Cold Reg. Sci. Technol., 20, 65–73, 1991. 

Vinkovic, I., Aguirre, C., Ayrault, M., and Simoëns, S.: Large-eddy Simulation of the Dispersion of Solid Particles in a Turbulent Boundary Layer, Bound.-Lay. Meteorol., 121, 283–311, 2006. 

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,, 2014. 

Xue, M., Droegemeier, K. K., Wong, V., Shapiro, A., Brewster, K., Carr, F., Weber, D., Liu, Y., and Wang, D.: The Advanced Regional Prediction System (ARPS) – A multi-scale nonhydrostatic atmospheric simulation and prediction tool. Part II: Model physics and applications, Meteorol. Atmos. Phys., 76, 143–165,, 2001. 

Zhang, J. and Huang, N.: Simulation of Snow Drift and the Effects of Snow Particles on Wind, Model. Simul. Eng., 2008, 408075,, 2008.  

Zwaaftink, C. D. G., 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, 2014. 

Short summary
Drifting snowstorms that are hundreds of meters in depth are reproduced using a large-eddy simulation model combined with a Lagrangian particle tracking method, which also exhibits obvious spatial structures following large-scale turbulent vortexes. The horizontal snow transport flux at high altitude, previously not observed, actually occupies a significant proportion of the total flux. Thus, previous models may largely underestimate the total mass flux and consequently snow sublimation.