the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

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

### Zhengshi Wang

### 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.

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.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 $\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}}={\left(\mathrm{\Delta}{x}_{\mathrm{1}}\mathrm{\Delta}{x}_{\mathrm{2}}\mathrm{\Delta}{x}_{\mathrm{3}}\right)}^{\mathrm{1}/\mathrm{3}}$,
in which Δ*x*_{i} is the grid spacing
along the streamwise (*i*=1), spanwise (*i*=2), and vertical directions
(*i*=3).
$\mathit{\rho}=p\left(\mathrm{1}-{q}_{\mathrm{v}}/\left(\mathit{\epsilon}+{q}_{\mathrm{v}}\right)\right)\left(\mathrm{1}+{q}_{\mathrm{v}}\right)/\left({R}_{\mathrm{d}}T\right)$
is the air density, in which *p*, *q*_{v}, *R*_{d}, 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. *u*_{i} is the
instantaneous wind speed component, and *x*_{i} is the position coordinate.
*t* is time, *δ*_{ij} is the Kronecker delta, $B=-g{\mathit{\rho}}^{\prime}/\mathit{\rho}$ is the buoyancy caused by the air density
perturbation *ρ*^{′}, and *g* is the acceleration due to gravity.
${p}^{*}={p}^{\prime}-\mathit{\alpha}\mathrm{\nabla}\left(\mathit{\rho}\mathit{u}\right)$ 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)

where ${\stackrel{\mathrm{\u0303}}{S}}_{ij}=\mathrm{0.5}\left(\partial {\stackrel{\mathrm{\u0303}}{u}}_{i}/\partial {x}_{j}+\partial {\stackrel{\mathrm{\u0303}}{u}}_{j}/\partial {x}_{i}\right)$ is the strain rate tensor and
$\left|\stackrel{\mathrm{\u0303}}{S}\right|=\sqrt{\mathrm{2}{\stackrel{\mathrm{\u0303}}{S}}_{ij}{\stackrel{\mathrm{\u0303}}{S}}_{ij}}$, *C*_{s} 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 10^{3} 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 ${x}_{{\mathrm{p}}_{i}}$ and ${u}_{{\mathrm{p}}_{i}}$ are the position coordinate and velocity of
the snow particle, respectively. *m*_{p} is the mass of the solid particle,
*V*_{r} is the relative speed between the snow particle and air, and
${T}_{\mathrm{p}}={\mathit{\rho}}_{\mathrm{p}}{d}_{\mathrm{p}}^{\mathrm{2}}/\mathrm{18}\mathit{\rho}\mathit{\nu}$ is the
particle relaxation time, where *ρ*_{p} is the particle density
(900 kg m^{−3}), *d*_{p} is the particle diameter and
$\mathit{\nu}=\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{5}}$ m^{2} s^{−1} is the kinematic viscosity of air.
*f*(*R**e*_{p}) can be expressed as (Clift et al., 1978)

where $R{e}_{\mathrm{p}}={V}_{\mathrm{r}}\mathrm{d}/\mathit{\nu}$ 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 ${V}_{{\mathrm{r}}_{i}}={\stackrel{\mathrm{\u0303}}{u}}_{i}\left({x}_{\mathrm{p}}\right)-{u}_{{\mathrm{p}}_{i}}+{u}_{i}^{\prime}$, in which ${\stackrel{\mathrm{\u0303}}{u}}_{i}\left({\mathit{x}}_{\mathrm{p}}\right)$ 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):

where ${T}_{\mathrm{L}}=\mathrm{4}\stackrel{\mathrm{\u0303}}{k}/\left(\mathrm{3}{C}_{\mathrm{0}}\stackrel{\mathrm{\u0303}}{\mathit{\epsilon}}\right)$
is the Lagrangian correlation timescale. Here, *C*_{0}=2.1 is the
Lagrangian constant, $\stackrel{\mathrm{\u0303}}{\mathit{\epsilon}}={C}_{\mathit{\epsilon}}{\stackrel{\mathrm{\u0303}}{k}}^{\mathrm{3}/\mathrm{2}}/\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}}$ 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 d*t*. $\stackrel{\mathrm{\u0303}}{k}$ is the subgrid turbulent kinetic energy and can
be obtained from the transport equation (Deardorff, 1980)

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

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)

in which *A*=0.2 is a constant, and ${\stackrel{\mathrm{\u203e}}{d}}_{\mathrm{p}}$ 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)

where *m*_{i} is the mass of particle and ${u}_{{\mathrm{p}}_{i\downarrow}}$ and
${u}_{{\mathrm{p}}_{i\uparrow}}$ 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

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

where *θ*_{i} and *v*_{i} are the impact velocity and angle,
respectively. Here, *θ*_{i} has a mean value of approximately
10^{∘} (Sugiura and Maeno, 2000), and 〈*v*_{i}〉 is
set to be the threshold of impact velocity. Considering
the steady-state saltation condition (one impact particle generates one
ejecta on average), 〈*v*_{i}〉 is determined by
setting the ejection number ${n}_{\mathrm{e}}=\mathrm{0.51}{v}_{\mathrm{i}}^{\mathrm{0.6}}{\mathit{\theta}}_{\mathrm{i}}^{\mathrm{0.16}}$ equal
to 1. In this way, the mean horizontal velocity of impact particles can be
obtained through $\u2329{u}_{{\mathrm{p}}_{i\downarrow}}\u232a=\u2329{v}_{\mathrm{i}}\u232a\mathrm{cos}\u2329{\mathit{\theta}}_{\mathrm{i}}\u232a$.

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):

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}_{*}(\equiv 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 10^{7} 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)

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).

## 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}_{*}=\mathrm{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).

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).

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.

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 *h*_{c} (referred as *Q*_{c}) to the total
suspension flux *Q*_{s} 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 *Q*_{c} to the total suspension flux is
non-negligible under various *h*_{c} values, the proportion of *Q*_{c} when
*h*_{c}=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 *Q*_{c} to the total suspension flux increases with friction velocity but
decreases with increasing *h*_{c}. From Fig. 8b, it can be seen that the
proportion of *Q*_{c} 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 $\left({\stackrel{\mathrm{\u0303}}{u}}_{\mathrm{i}}\left({\mathit{x}}_{\mathrm{p}}\right)\right)$ 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.

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.

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.

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.

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.

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 (http://www.caps.ou.edu/ARPS/, Xue et al., 2001).

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.

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, https://doi.org/10.1103/PhysRevLett.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, https://doi.org/10.5194/tc-11-3011-2017, 2017.

Huang, N. and Wang, Z.: A 3-D simulation of drifting snow in the turbulent boundary layer, The Cryosphere Discuss., 9, 301–331, https://doi.org/10.5194/tcd-9-301-2015, 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, https://doi.org/10.5194/acp-16-7523-2016, 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, https://doi.org/10.1029/2010JD015419, 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, https://doi.org/10.1029/2002jd003327, 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, https://doi.org/10.1029/2004jd004657, 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, https://doi.org/10.1029/2011JD015828, 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, https://doi.org/10.5194/tc-8-395-2014, 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, https://doi.org/10.1007/s007030170027, 2001.

Zhang, J. and Huang, N.: Simulation of Snow Drift and the Effects of Snow Particles on Wind, Model. Simul. Eng., 2008, 408075, https://doi.org/10.1155/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.