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

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.


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 Published by Copernicus Publications on behalf of the European Geosciences Union.
Z. Wang and S. Jia: A simulation of a large-scale drifting snowstorm in the turbulent boundary layer 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 andWang, 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 largeeddy 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 smallscale 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.

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 nonhydrostatic Navier-Stokes equation is solved (Xue et al., 2001): where "∼" represents variables that are filtered and the filtering scale is = ( x 1 x 2 x 3 ) 1/3 , in which x i is the grid spacing along the streamwise (i = 1), spanwise (i = 2), and vertical directions is the air density, in which p, q v , R d , and T are the pressure, specific humidity, gas constant (287.0J 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ρ /ρ 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) where S ij = 0.5 ∂ u i /∂x j + ∂ u j /∂x i is the strain rate tensor and | S| = 2 S ij S ij , C s is the Smagorinsky coefficient that is determined locally by the dynamic Lagrangian model (Meneveau et al., 1996).

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 p i and u 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 p = ρ p d 2 p /18ρν is the particle relaxation time, where ρ p is the particle density (900 kg m −3 ), d p is the particle diameter and ν = 1.5 × 10 −5 m 2 s −1 is the kinematic viscosity of air.f Re p can be expressed as (Clift et al., 1978) where Re p = V r d/ν 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 r i = u i x p − u p i + u i , in which u i x p 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 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) where θ is the potential temperature and θ 0 is the surface potential temperature.

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 particleborne 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 d 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 p i↓ and u p i↑ 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 liftoff 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 e = 0.51v 0.6 i θ 0.16 i equal to 1.In this way, the mean horizontal velocity of impact particles can be obtained through 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.

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 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 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 Results and discussions

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).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 mea-  surements (Nishimura and Nemoto, 2005;Schmidt, 1982).Therefore, the established model is able to produce a largescale 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.

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 al- most 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, 1982Schmidt, , 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, 1982Schmidt, , 1984;;Tabler, 1991).

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 u i x p 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 smallscale 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    Figure 10.Evolution 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).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.

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

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

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

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

Figure 4 .
Figure 4. Variation in average particle diameter versus height.

Figure 5 .
Figure 5. Particle size distribution at various heights.

Figure 6 .
Figure 6.Variations in snow transport flux versus height.

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

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

Figure 11 .
Figure 11.Time evolutions of the thickness of the drifting snowstorm under various friction velocities.
www.the-cryosphere.net/12/3841/2018/The Cryosphere, 12, 3841-3851, 2018 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.