Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

Contents lists available at ScienceDirect

Journal of Wind Engineering & Industrial Aerodynamics
journal homepage: www.elsevier.com/locate/jweia

Structure of drifting snow simulated by Lagrangian particle dispersion
model coupled with large-eddy simulation using the lattice
Boltzmann method
Tsutomu Watanabe a, *, Shuhei Ishikawa b, Masayuki Kawashima a, Kou Shimoyama a,
Naoyuki Onodera c, Yuta Hasegawa c, Atsushi Inagaki daInstitute of Low Temperature Science, Hokkaido University, Sapporo, 060-0819, Japan
Graduate School of Environmental Science, Hokkaido University, Sapporo, 060-0810, Japanc
Center for Computational Science & e-Systems, Japan Atomic Energy Agency, Kashiwa, 277-0871, Japand
Department of Transdisciplinary Science and Engineering, Tokyo Institute of Technology, Tokyo, 152-8552, Japanb

A R T I C L E I N F O

A B S T R A C T

Keywords:
Saltation layer height
Force balance
Transition probability
Snow streamers
Low-speed streaks
Ejection
Sweep
Coherent turbulent motion
Atmospheric surface layerThis paper describes simulations of drifting snow in the atmospheric surface layer. We develop a Lagrangianparticle dispersion model coupled with a large-eddy simulation code based on the central-moment lattice
Boltzmann method. The model reproduces typical features of drifting snow observed in the field, such as thedependency of the mass transport rate on the flow velocity, the kink in the vertical mass flux profile near thesaltation layer height, and the variations in particle size distribution with the flow velocity and height. The
saltation layer height determined directly from the net forces acting on the airborne particles is found to increase
monotonically with increasing flow velocity, unlike conventional estimates, which tend to saturate as the flow
velocity increases. Using the vertical transition probabilities of individual particles, the transition from saltation
to suspension is confirmed to occur near the estimated saltation layer height. A composite analysis shows thatsnow streamers (dense particle clouds elongated in the streamwise direction and meandering laterally in thesaltation layers) are closely associated with small-scale low-speed streaks in the near-surface flows. Particularly
dense snow streamers are more likely to occur around streaks modulated by high-speed coherent flows of muchlarger spatial scales.

1. Introduction
Drifting snow is a significant factor in hazard control because it re­duces visibility, creates obstructions through snowdrifts, and producessnow cornices that may trigger avalanches. The large-scale redistribu­
tion of snow is known to affect the surface water balance in polar regions
and hilly terrain. Terrestrial ecosystems and underground frozen soil inhigh-latitude and/or high-altitude areas, which are highly sensitive tothe spatial distribution of vegetation and micro-topography, can also be
profoundly affected by the micro-scale environments produced by het­
erogeneous depositions of snow (see Mott et al., 2018 for a compre­
hensive review of wind-driven snow cover dynamics).
Drifting snow is essentially a non-stationary two-phase flow, inwhich airborne particles and turbulence strongly interact. Further
complexity is added by the spatio-temporal intermittencies induced bycoherent turbulent motion onto the near-surface flow velocities and thesurface processes that feed particles into the flow (Doorschot et al.,
2004; Aksamit and Pomeroy, 2018a, 2018b). Therefore, besides thepivotal role of observational and experimental studies, the need fornumerical modelling cannot be overemphasized as a means of compre­
hending the nature of drifting snow phenomena that occur in the at­mospheric surface layers.
The first attempt at three-dimensional numerical modelling was
made by Uematsu et al. (1991). They modelled drifting snow using an
Eulerian advection–diffusion equation for the bulk density of driftingsnow particles, which were driven by an airflow simulated by the
Reynolds-averaged Navier–Stokes (RANS) equations. Their model wasapplied to the simulation of snowdrifts around obstacles (Sato et al.,1993). Since then, Tominaga and Mochida (1999), Beyers et al. (2004),
and many subsequent studies have refined this type of model and* Corresponding author.
E-mail address: t-wata@lowtem.hokudai.ac.jp (T. Watanabe).
https://doi.org/10.1016/j.jweia.2024.105783
Received 15 November 2023; Received in revised form 23 May 2024; Accepted 23 May 2024
Available online 1 June 2024
0167-6105/© 2024 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

developed improved versions (see Tominaga, 2018; Zhou and Zhang,2023 for comprehensive reviews of numerical studies on snowdrifts
around obstacles). Owing to their computational efficiency, Eulerianmodels are highly suitable for regional-scale studies (Liston and Sturm,
1998; Gauer, 2001; Liston et al., 2007; Tominaga et al., 2011) and for
operational purposes (Pomeroy et al., 1993).
A Lagrangian approach, in which the trajectories of individual par­ticles are traced by integrating the equation of motion, provides a more
explicit representation of the particles’ behaviour in the flow. The studyof Shao and Li (1999) was probably the first to use a Lagrangian particle
model coupled with a RANS airflow model for simulating drifting snow.
Nemoto and Nishimura (2004) took a similar approach in their simu­
lations of two-dimensional snow-particle trajectories in a coupled
one-dimensional flow, and this model was later used by Niiya andNishimura (2017, 2022). Recently, this type of model has been applied
to three-dimensional simulations of drifting snow over roofs (Chen et al.,2021; Liu et al., 2022; Zhou et al., 2023). A remarkable feature of the
Lagrangian approach is its ability to realistically reproduce the parti­cle–surface interactions. This is because the key surface processes of the
aerodynamic entrainment of particles, the rebound of impacting parti­cles, and the splash entrainment of the snow-bed particles can bedirectly incorporated into the tracing of individual trajectories.
The intermittency induced by turbulent coherent motion, which
RANS models are unable to represent, has motivated the development of
next-generation models, namely Lagrangian particle models coupledwith large-eddy simulation (LES) models of atmospheric surface-layer
flow. Pioneering studies were conducted by Dupont et al. (2013) for
drifting soil particles and by Groot Zwaaftink et al. (2014) for driftingsnow. Recent studies have used LES-coupled Lagrangian snow-particle
models to investigate the horizontal structure of drifting snow (Huang
and Wang, 2016), the intermittent nature of drifting snow (Okaze et al.,2018), large-scale drifting snowstorms in the atmospheric boundary
layer (Wang and Jia, 2018), sublimation of saltating snow particles(Sharma et al., 2018; Wang et al., 2019), heat and water-vapour transfer
during drifting snow events (Sigmund et al., 2022), and the dependencyof drifting snow on the snow properties (Melo et al., 2022).
Besides solving the Navier–Stokes equations directly, the lattice
Boltzmann method (or algorithm) offers another approach to flow
simulation (Chen and Doolen, 1998; Aidun and Clausen, 2010; Krüger
et al., 2017). This method considers a fluid as an ensemble of micro­
scopic fluid particles moving around and colliding with each other, anddescribes fluid motion by the spatio-temporal evolution of the velocity
distribution function of the fluid particles. The macroscopic fluid den­sity, velocity, and pressure derived from low-order moments of the
distribution function satisfy the Navier–Stokes equations. The method,which has been used successfully in applications ranging fromisothermal single-phase flows to fully compressible multiphase flows,
can be considered an established solver for the Navier–Stokes equations.
As noted by Kareem (2020), it is gaining popularity in engineering ap­
plications, such as simulating air flow around buildings (Lenz et al.,2019; Han et al., 2020) and ships (Syms, 2008), and aerodynamics ofcomplex structures (Andre et al., 2015; Islam et al., 2017; Buffa et al.,2021). Recently, the lattice Boltzmann method was confirmed toperform well in simulating atmospheric boundary-layer flows (Ahmad
et al., 2017; Inagaki et al., 2017; Onodera et al., 2021; Feng et al., 2021)
and atmospheric surface-layer flows including plant canopies (Wata­nabe et al., 2020, 2021). The computational parallelism shared by thelattice Boltzmann algorithm and the Lagrangian particle model makes
their combination appropriate for high-performance computations ofdrifting snow. The algorithm’s flexibility in terms of incorporatingcomplex boundaries is advantageous for simulations of drifting snow
over obstacles with complex geometries or time-varying micro-­
topographies (e.g., snowdrifts) produced by erosion and deposition, andthus will be beneficial for future applications for snow engineering. The
application of the lattice Boltzmann method to the simulation of driftingsnow was initiated by Tanji et al. (2021), who focused on simulating thewind field around various type of fences and the resulting snowdrift
formation. For that purpose, their model did not consider particle–sur­face interactions or the effects of particle drag on the flow, making it
unable to reproduce the self-sustaining nature of drifting snow.
In this paper, we first describe a new Lagrangian particle dispersion
model that is fully coupled with the lattice Boltzmann method and in­cludes particle–surface interaction processes. The adopted lattice
Boltzmann code, which uses the central-moment-based multi-relaxa­
tion-time collision scheme (Geier et al., 2006, 2015), is the same as thatused by Watanabe et al. (2020, 2021) to simulate canopy turbulence,
and achieves performance comparable to a Navier–Stokes model(Watanabe, 2004, 2009). Aerodynamic entrainment, rebounding, andsplash entrainment are included as particle-feeding processes on thesnow bed, similar to previous drifting-snow models. We examine themodel’s ability to reproduce the properties of drifting snow, including
instantaneous and statistical features, and investigate the dynamics ofsaltating particles. Finally, we elucidate the relationship between thehorizontal structures of drifting snow and turbulence in the vicinity ofthe surface.
2. Model description
A Cartesian coordinate system is adopted, in which x, y, and z denotethe streamwise, lateral, and vertical directions, respectively. Both thevector notation and the component form of a vector, i.e., x = (x,y,z), areused in this paper.
2.1. Large-eddy simulation of the near-surface flow
The model used for simulating the turbulent flow in the atmosphericsurface layer has been described in our previous reports (Watanabe
et al., 2020, 2021); a summary is given here. The model is based on thecentral-moment lattice Boltzmann method in a three-dimensional 27-ve­
locity (D3Q27) lattice system with a uniform lattice spacing Δx. The
method predicts the spatio-temporal evolution of the distributionfunction fijk (x, t) of the microscopic fluid-particle velocities (ic, jc, kc),
where the indices i, j, k ∈ {− 1, 0, 1} represent the directions of dis­cretized velocities and c = Δx/Δt denotes the lattice velocity unit, withΔt being the timestep. The velocity distribution function is updated ateach timestep through a collision–forcing step followed by a streaming
step.
The collision–forcing step begins with the central moments of the
distribution function being evaluated byκpqr (x, t) =

1 ∑
1
1
∑
∑
i=− 1 j=− 1 k=− 1(ic − u)p (jc − v)q (kc − w)r fijk (x, t),(1)where κpqr denotes the central moment of the directional order p, q, r ∈{0, 1, 2}, and u, v, and w are the streamwise, lateral, and vertical com­ponents of the macroscopic velocity u = (u,v,w), respectively. Using themulti-relaxation-time collision scheme and a forcing scheme applied in
the central-moment space, the post-collisional central moments κ*pqr are
calculated in terms of the pre-collisional central moments, the macro­scopic fluid density ρ, the macroscopic velocity, a net volumetric force()
F = Fx , Fy , Fz , and the collisional relaxation coefficients ω1 , ⋯, ω10 asκ*pqr (x, t) = func(κ, ρ, u, FΔt, ω1 , ⋯, ω10 ).

(2)

Details of the above relationship are described in Watanabe et al.
(2020). The inverse transform of Eq. (1) is then applied to κ*pqr to*calculate the post-collisional distribution function fijk
, which is trans­
lated to the adjacent grid points in the next streaming step as*fijk (x + icΔt, y + jcΔt, z + kcΔt, t + Δt) = fijk(x, y, z, t).

(3)

After appropriate boundary conditions have been applied to this new2

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

distribution function, the macroscopic variables at the next timestep are
evaluated from the low-order raw moments as

ρ=

1 ∑
1
1
∑
∑
i=− 1 j=− 1 k=− 1

ρu =

ρv =

fijk ,1 ∑1
1∑∑
i=− 1 j=− 1 k=− 1

(7a)

(4a)

dxp
= up ,
dt(7b)

where ρp is the particle density, |Δu| is the magnitude of the velocity
difference between the particle and the flow, us is the subgrid-scale flow
velocity experienced by the particle, g is the acceleration due to gravity,̂ is the vertical unit vector. The particle’s drag coefficient cd isand k(ic)fijk +

Fx Δt
,
2

(4b)

(jc)fijk +

Fy Δt
,
2

(4c)

1 ∑
1
1
∑
∑
i=− 1 j=− 1 k=− 1)(
dup 3 cd ρ̂=|Δu| u + us − up − g k,
dt 4 dp ρp

1 ∑
1
1
∑
∑

Fz Δt
,
ρw =
(kc)fijk +2
i=− 1 j=− 1 k=− 1

ρ

p= ,
3

evaluated by the empirical formula (White, 1974; Okaze et al., 2018)
cd =

(4d)

(8)

where Rep = dp |Δu|/ν is the particle Reynolds number. This expression
is valid for the range 0 ≤ Rep ≤ 2 × 105 (White, 1974). The flow velocity
at the particle position is calculated from the simulated velocities atsurrounding grid points by interpolating bilinearly in the horizontal
directions and logarithmically in the vertical direction. For the particles
̃ (zb ) =below the lowest grid level zb , the flow velocity is given by u(z) /ũ (zb ) is the horizontalln[(z + z0s ) /z0s ] /ln[(zb + z0s ) /z0s ], where ubilinear interpolation of the flow velocity at z = zb . The subgrid-scale
flow velocity us is estimated based on the Lagrangian stochastic model
described by Weil et al. (2004), as detailed in Appendix A. It should be
noted that while these treatments (interpolation of the resolved veloc­
ities and modelling of the unresolved velocities) attempt to alleviate theproblem caused by insufficient near-surface resolution, some of the re­sults presented in subsequent sections may be dependent on thesetreatments. However, as the Stokes number near the surface is quite
large (see Table 1), the effect of the flow uncertainty on the near-surface
particle trajectories may not be significant. Comparisons with muchhigher resolution simulations are needed to confirm this and are left for
future research.
The first term on the right-hand side of Eq. (7a) (denoted here as ap )represents the acceleration induced by the aerodynamic force, the re­action of which exerts a drag force on the flow. Hence, Fp in the flow
model described above is given by(4e)where p is the static pressure. To enable simulations of high Reynolds
number flows, the subgrid eddy viscosity νe is parametrized using thecoherent-structure Smagorinsky model of Kobayashi (2005) and incor­
porated into the collisional relaxation coefficients for the even-order
central moments as 1/ω = 3(ν + νe ) + 1/2, where ν is the molecular
kinematic viscosity (Watanabe et al., 2020).
The net volumetric force F is given by the sum of the flow-driving
force Fdrv , the surface friction Fs , and the particle drag Fp . The driving
force of the flow is given by)(ρu2*f ̂Fdrv =
i,(5)
ztop
where u*f is the friction velocity before the onset of drifting (hereafterreferred to as the freestream friction velocity), which is specified foreach simulation run, ̂i is the streamwise unit vector, and ztop is theheight of the computational domain. The surface friction is modelled bya modified bulk drag equation that accounts for the momentum sinkcaused by particles drifting below the lowest grid level as⃒⃒/
1
Fs = − ρCDs ⃒uh,b ⃒uh,b Δx − Fp,b ,224
6√̅̅̅̅̅̅̅̅ + 0.4,+
Rep 1 + Rep(6)

Fp = −

where CDs is the bulk drag coefficient of the snow surface, uh is thehorizontal component of the flow velocity, and the subscript b repre­
sents variables at the lowest grid level. The bulk drag coefficient isevaluated assuming a logarithmic velocity profile with a prescribedroughness length z0s . The particle drag force Fp is the sum of the dragforces exerted by all particles within each grid cell, as described in thenext subsection.

∑

ρp

πd3p
6

(9)

ap ,

where the summation is taken over all particles in each grid cell.
2.2.2. Aerodynamic entrainment
The entrainment rate Nae [m− 2 s− 1], defined as the number of par­ticles picked up by the flow per unit time and unit area, follows theexcess shear stress rule (Anderson and Haff, 1991)()(10)
Nae = ηρ u2* − u2*t .

2.2. Lagrangian model of snow particle motion
The drifting snow particles are assumed to be spherical, and theirvolume and frontal area are calculated from the particle diameter dp .
Drifting snow particles, while initially complex in shape, tend to breakinto crystal fragments over the course of multiple impacts on the bed(Sato et al., 2008; Comola et al., 2017), and many of which naturally
become more or less rounded rather than retain their original shape(although perfect spheres are rarely observed) (Gordon and Taylor,2009). The sublimation of the particles is not considered.

Table 1
Relevant parameters for each case of simulations. u1m is the simulated mean
wind speed at a height of 1 m. The Reynolds numbers for the flow and particles
are defined as Reτ = u*f L/υ and Rep,max = u1m 〈dp 〉/υ, respectively, where L (=4m) is the height of the domain. The bulk Stokes number St and the viscous Stokes
number

St +
(

are
)

definedas

St = ρp 〈dp 〉2 u1m /(18ρνL)ρp 〈dp 〉2 u*f 2 / 18ρν2 , respectively.

2.2.1. Equation of motion
The particles, once seeded to the flow, are accelerated by aero­dynamic drag and gravity, and so their velocities and positions changewith time. The particle velocity up and the particle position xp aredetermined by
3

andSt + =

u*f [m s− 1]

u1m [m s− 1]Reτ

Rep,max

St

St+

0.3
0.35
0.4
0.5
0.6
0.8

6.3
7.2
8.0
9.3
10.5
12.7

9.6 × 104
1.1 × 105
1.3 × 105
1.6 × 105
1.9 × 105
2.6 × 105

102
116
129
150
169
205

0.19
0.22
0.26
0.29
0.32
0.39

8.9 × 102
1.2 × 103
1.6 × 103
2.5 × 103
3.6 × 103
6.3 × 103

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
Here, the friction velocity is calculated from u* = |Fs |Δx/ρ, and thethreshold friction velocity u*t is estimated by (Bagnold, 1941)√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
ρp − ρ
u*t = A(11)
gdp ,and momentum as (Gauer, 2001; Comola and Lehning, 2017)
Nspl = min(Nse , Nsm ),[(ρη = kae4 ρs,
π〈dp 〉2 ρp

⃒
⃒2
)
2
1 − δrbd rrbd
− ϵf mimp ⃒uimp ⃒
,0 ,
⃒ ⃒2
〈mspl 〉〈⃒uspl ⃒ 〉 + 2ϕ

(15b)

)
⎤
⎡ (
1 − δrbd rrbd − μf mimp uimp⎣⃒⃒
Nsm = max
, 0⎦,
〈mspl 〉〈⃒uspl ⃒〉〈cos αspl 〉〈sin βspl 〉(15c)

Nse = maxwith the constant parameter A set to 0.18 (Clifton et al., 2006). The
entrainment coefficient η [N− 1 s− 1] is given by (Clifton and Lehning,
2008)
(12)

(15a)
]

where Nse and Nsm are the numbers estimated from energy and mo­
mentum conservation, respectively. Here, δrbd = 1 if rebound occurs andδrbd = 0 otherwise, ϵf (=0.66; Ammi et al., 2009) and μf (=0.35; Riceet al., 1995) are the fractions of kinetic energy and momentum dissi­
pated during the impact, respectively, mimp is the mass of the impacted
particle, mspl , uspl , αspl , βspl are the mass, velocity, and vertical andhorizontal angles of the splashed particle, respectively, and ϕ [J] is theinter-particle cohesion energy, as determined below. The product 〈cosαspl 〉〈sinβspl 〉 is set to be constant at 0.8 (Comola and Lehning, 2017).
The initial motion of a splashed particle is defined by the initial speed,
⃒
⃒
⃒
⃒0.3
which follows an exponential distribution with 〈⃒uspl ⃒〉 = 0.25⃒uimp ⃒where kae [m2 N− 1 s− 1] is a constant parameter set to 1, ρs is the surface
bulk density of snowpack (set to 150 kg m− 3), and angle brackets denotethe mean value. The size of each entrained particle is randomly chosenfrom the size distribution of the snow-bed particles, as specified later.
Note that ρs is included in Eq. (12) to account for the number of particles
available on a unit area of the snow-bed surface and should inherently
depend on the snow-bed conditions. However, in most cases, theentrained particles make up only a small fraction of drifting particles(see Section 3.2), so the results presented in this paper are not sensitive
to this value.
The initial speed vae and vertical angle αae of the entrained particle
are specified following Clifton and Lehning (2008) in a similar manner
to Groot Zwaaftink et al. (2014). The magnitude of the initial velocity ofthe entrained particles follows a log-normal distribution with 〈vae 〉 =3.3u* and standard deviation 0.5〈vae 〉. The mean angle 〈αae 〉 depends on
the particle size and is given by
)]}
{
[
(
π
dp
〈αae 〉 =,(13)
75 − 55 1 − exp −180
1.75 × 10− 4(Anderson and Haff, 1991), the vertical angle, which obeys an expo­
nential distribution with 〈αspl 〉 = 50◦ (Kok and Renno, 2009), and thehorizontal angle, which is chosen from a normal distribution around theimpact direction with a standard deviation of 15◦ (Xing and He, 2013).
The size of the splashed particle is randomly taken from the size dis­tribution of the snow-bed particles.
Fig. 1 shows the relationship between the impact speed and thenumber of ejected particles (including both rebounded and splashedparticles) obtained from Monte Carlo simulations of the particles’
impact–rebound and splashing for different values of ϕ. We followed theprocedure of Comola and Lehning (2017), except in using the reboundand splash models described above and the particle size distribution
defined by Eq. (16). As the results for ϕ = 3 × 10− 9 J best fit thewind-tunnel measurements reported by Sugiura and Maeno (2000), thisis taken as the standard value. However, considering the uncertainty inϕ, which varies depending on the snow-bed conditions, non-cohesive (ϕwhere individual realizations of the angle obey a log-normal distri­bution with a standard deviation of 0.5〈αae 〉. The horizontal angle ofejection is assumed to be parallel to the flow. The entrained particles arelaunched from a height z = (Δt/2)vae sinαae at a randomly chosen hori­zontal position within each grid area.
2.2.3. ReboundSnow particles that impact the surface may rebound with a proba­
bility Prbd defined by (Anderson and Haff, 1991; Groot Zwaaftink et al.,
2014)
⃒)]
⃒
[
(
Prbd = 0.9 1 − exp − 2.0⃒uimp ⃒ ,(14)where uimp [m s− 1] is the impact velocity. The restitution coefficient of
the rebound rrbd is assumed to obey a normal distribution with mean
0.55 and standard deviation 0.1, based on the measurements of Araoka
and Maeno (1981), Kosugi et al. (1995), and Nishimura and Hunt⃒⃒
(2000). The rebound speed is calculated as rrbd ⃒uimp ⃒, and the verticalangle of rebound follows an exponential distribution with a mean valueof 45◦ (Kok and Renno, 2009), while the horizontal angle is parallel tothe incident particle.
2.2.4. Splash
The splash entrainment is the most unknown of the particle feeding
processes. The number of splashed particles can be estimated relatively
reasonably well under the constraints of conservation of kinetic energy
and momentum. However, there are only a few examples of systematic
measurements of particle splashing speeds and angles on snow surfaces(e.g., Sugiura and Maeno, 2000), even though they should be affected bymany factors, such as snow-bed conditions, properties of impacted
particles, impact speed, impact angle. Therefore, statistical relationships
established for soil and sand particles from many laboratory and nu­
merical studies are substituted here.
The number of particles that are splashed from the snow bed upon asingle particle impact Nspl is modelled based on conservation of energy

Fig. 1. Mean number of particles ejected (rebounded and splashed) upon asingle impact calculated by Monte Carlo simulations using Eqs. (14)–(16) withdifferent cohesive energies. Wind-tunnel measurements reported by Sugiuraand Maeno (2000) are also plotted for comparison, with error bars indicating
standard deviations. Each simulation performed at 0.1 m s− 1 intervals in impact
speed involves 106 particles impacting the surface at randomly sampled angles
ranging from 5 to 15◦ . The sizes of impacting and splashing particles are
randomly chosen from the gamma distribution defined by Eq. (16).
4T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783= 0) and highly cohesive (ϕ = 1 × 10− 8 J) conditions are also consideredin the following simulations. Incidentally, when the same Monte Carlosimulation is performed for single-size non-cohesive particles (e.g., ho­⃒⃒
mogeneous sand), the relationship Nspl ≈ 0.5⃒uimp ⃒ is observed. This isconsistent with previous measurements and numerical simulations (seeFig. 3a of Comola and Lehning, 2017).
2.2.5. Numerical implementation
Aerodynamic entrainment is invoked at a constant time interval Δtae ,with the number of entrained particles calculated from Eq. (10) using
the shear stress ρu2* at each surface grid point averaged over Δtae . By
considering the impact–rebound process every timestep, the kinetic
energy and momentum remaining after all impact events that occur ineach grid area are accumulated over the time interval Δtspl , and theseaccumulated amounts are used to calculate the number of splashedparticles at each Δtspl . Particles entrained (or splashed) at the same timeand location are treated as a group of Nae (or Nspl ) single-sized particles
that share the same velocity and position, and each group’s trajectory is
tracked thereafter. The number of particles belonging to each group (Naeand Nspl ) retains the fractional part of the number calculated from Eqs.
(10) and (15), which is ignored if the number is less than 1. The model
code is written in the NVIDIA CUDA language, and the individual tra­jectories are computed in parallel using graphics processing units.

Fig. 3. Horizontal distribution of particle groups obtained at the same instantas for Fig. 2 for the case of u*f = 0.5 m s− 1, shown over the half-width of thecomputational domain.

time. The timestep is the same for the flow and particle motions, anddecreases linearly with the freestream friction velocity: from 0.236 ms
for u*f = 0.3 m s− 1 to 8.84 × 10− 2 ms for u*f = 0.8 m s− 1. The timeintervals of aerodynamic entrainment and splash entrainment are set toΔtae = Δtspl = 86Δt, which corresponds to 20 ms for the case of u*f =0.3 m s− 1. With these settings, the number of particle groups tracked
ranges from 106 to 107 depending on u*f , with each group typicallycontaining 100 to 103 particles. During the simulations, snow particles
moving upward beyond the top boundary and stopping (not rebound­ing) at the surface are discarded.
The density and kinematic viscosity of the air are set to ρ = 1.34 kg
m− 3 and ν = 1.24 × 10− 5 m2 s− 1, respectively, and the density of snow
particles is set to ρp = 918 kg m− 3. Following Nemoto and Nishimura(2004), the size distribution of the erodible particles on the snow bed isgiven by a gamma distribution with 〈dp 〉 = 2 × 10− 4 m, expressed as()( )dαp− 1dp
f dp = α
exp −
,
(16)
β Γ(α)
β2.3. Simulation setup
A computational domain of 12.8 × 12.8 × 4.0 m3 with a uniform grid
resolution of 0.05 m is used to simulate drifting snow over a flat andhomogeneous snow surface under six different wind conditions of u*f =0.3, 0.35, 0.4, 0.5, 0.6, and 0.8 m s− 1. The top boundary at ztop = 4 m is arigid free-slip wall for the airflow, while the bottom boundary (snowsurface) is a rigid wall with friction, as explained in Section 2.1. Theserigid-wall boundaries are represented in the lattice Boltzmann algorithm
by the half-way specular reflection of the wall-ward components of thevelocity distribution function (Watanabe et al., 2020, 2021). The hori­zontal boundaries are set to be periodic to mimic a developed
surface-layer flow over a homogeneous surface. The flow simulations
are initialized with a uniform streamwise velocity field withsmall-magnitude random perturbations of the vertical velocity compo­
nent. After the turbulence statistics (e.g., mean flow velocity, Reynolds
stress, turbulent kinetic energy) reach a steady state, the snow-particle
processes are activated to permit the drifting snow to develop overwhere α = 5 and β = 4 × 10− 5 m are constant parameters and Γ(α)
denotes the gamma function. The mean particle mass used in Eq. (15)can be evaluated using 〈d3p 〉 = α(α + 1)(α + 2)β3 . Prior to the simula­tions, a look-up table of cumulative probabilities for 2 × 10− 5 m in­tervals up to the maximum diameter of 6.2 × 10− 4 m is prescribed using
this formula. The roughness length of the surface is assumed to remain
constant at z0s = 2.3 × 10− 4 m (Kondo and Yamazawa, 1986; Gromke
et al., 2011). Table 1 summarizes the parameters related to the flow andparticles simulated with these settings. Obviously, the flow is fully tur­bulent, with conditions such that the particles behave as tracers in mostof the domain (if the gravity is negligible) and as ballistic projectiles in
the region very close to the surface.
Preliminary simulations confirmed that the total mass transport rate
reaches steady state by 200 s from the onset of the drifting. Therefore,
datasets for post-processing are saved every 1 s for 200–1000 s. Aver­
aging these 800 datasets has also been confirmed sufficient to obtainstable statistics on the fully developed state in which the flow–particle
interactions are in equilibrium at all heights in the domain. Unlessotherwise stated, the mean values at a given height are calculated astime averages of the horizontal mean values obtained from all datasets,
and are denoted by an overbar. Additionally, the particle-related sta­tistics are evaluated at logarithmically equidistant heights, with a finer
resolution near the surface than the flow statistics evaluated at each gridlevel.

Fig. 2. Vertical–streamwise distribution of particle groups 10 s after the onsetof drifting for different values of u*f : (a) 0.35, (b) 0.5, and (c) 0.8 m s− 1. All
particle groups (see Section 2.2.5) in the lower half of the computational
domain are rendered as dots with opacity linearly related to the total crosssectional area of the particles in each group.
5

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

3. Results and discussion
3.1. Instantaneous visualization
For a visual inspection of the simulated drifting snow, Figs. 2 and 3
present the vertical and horizontal distributions of the particle groups
obtained 10 s after the onset of drifting. Fig. 2 (and its animated version,
provided as Supplementary Video 1) illustrates the differences in theheight and density of drifting snow (and the resulting visibility reduc­tion) and their growth rates over time at different flow velocities. Fig. 2
also displays the spatial density patterns tilting downstream in the nearsurface region, which are affected by coherent eddies developing in thesheared surface-layer flow. The horizontally fluctuating structurescommonly seen in the aeolian particle transfer phenomena occurring
near the surface (Bass, 2008) are illustrated in Fig. 3 (their temporalvariations can be seen in the animation provided as Supplementary
Video 2). The spatial relationship between the horizontal patterns in theconcentration of drifting particles and the instantaneous flow velocities
near the surface will be discussed in Section 3.3.
3.2. Statistics and vertical structure of the drifting snow
3.2.1. Total mass transport rate
The total mass transport rate Q [kg m− 1 s− 1], calculated by inte­
grating the streamwise mass flux with respect to height, is comparedwith wind-tunnel measurements and empirical formulas derived fromsuch measurements in Fig. 4a. All these wind-tunnel experiments wereconducted on smoothly conditioned snow surfaces. Here, for Sørensen’s(2004) formula, the parameters proposed by Vionnet et al. (2014) areused and the friction velocity at the impact threshold is assumed to be
0.26 m s− 1, which is the equilibrium surface friction velocity obtained
from the simulation with ϕ = 0 for the lowest u*f (=0.3 m s− 1).
The simulated values of the transport rate at equilibrium shown in
Fig. 4a (filled symbols) are generally smaller than those measured inwind tunnels, especially for small u*f cases; this trend is similar to thatreported by Melo et al. (2022). This trend at low u*f may be partly
attributable to the intermittency of the aerodynamic entrainment. In thesimulations for low-wind-speed conditions, the surface friction veloc­
ities do not exceed the threshold simultaneously at all locations becauseof the higher thresholds for larger particles and the heterogeneity ofnear-surface turbulent flow. In such a situation, drifting snow occurs
only where coherent eddies from aloft enhance the aerodynamic
entrainment followed by the splash entrainment of snow-bed particles,
as reported by Okaze et al. (2018) and as will be discussed in Section 3.3.
This situation contrasts with the wind-tunnel experiments conductedunder relatively uniform wind conditions. Meanwhile, even in the
non-cohesive case (ϕ = 0), where the modelled mass transport rate ishighest, halving the entrainment threshold (Eq. (11)) or doubling theentrainment coefficient (Eq. (12)) only increases the transport rate by
35% and 20% at u*f = 0.3 m s− 1, respectively, and does not compensate
for the difference from the wind-tunnel results (figure not shown). Thisis likely due to the negative feedback between wind and particles; as thenumber of airborne particles increases, wind speed near the surface
decreases and the entrainment is reduced.
More importantly, the large sensitivity of Q to ϕ in the low-u*f region
indicates the importance of accurate parameterization for the splash
entrainment at low wind speeds. As shown in Fig. 1, ϕ = 3 × 10− 9 J is thevalue optimized through comparison with direct measurements of the
impact–splash process on a snow bed (Sugiura and Maeno, 2000), but
the model results using this value nevertheless underestimate thetransport rates measured at low wind speeds. The wind-tunnel experi­
ment of Sugiura and Maeno (2000) was conducted at − 15 ◦ C to avoid
rapid sintering between snow particles, but the actual cohesivenessduring the experiment is unknown. Indeed, the model results for ϕ =0 are close to Sørensen’s formula, but the impact threshold used for thisformula is higher than that (0.2 m s− 1) measured by Nishimura and HuntFig. 4. Total mass transport rates simulated using the different cohesive en­ergies: (a) comparison with wind-tunnel measurements plotted versus thefreestream friction velocities; (b) comparison with field observations plotted asa function of the mean wind speed at a height of 1 m. In (a), some model results
have been shifted horizontally for clarity; the edges of the error bars indicate
the initial transport rates averaged for the first 30 s of each simulation; anasterisk denotes the initial transport rate calculated using three times the Nspl
value for ϕ = 0; open circles and open triangles denote wind-tunnel measure­ments by Sugiura et al. (1998) and Okaze et al. (2012), respectively; the lines
represent the empirical formulas of Nishimura and Hunt (2000) and Sørensen(2004) for snow, and Bagnold’s (1941) formula for sand; grey shading indicatesthe range of u*t for the snow-bed particles, and the vertical long-dashed line
indicates the u*t corresponding to the mean diameter. In (b), asterisks denotethe transport rates observed by Kobayashi et al. (1970) and filled circles andopen triangles represent those observed in the saturation zone and in theneutral zone, respectively, by Takeuchi et al. (1975) and Takeuchi (1980).

(2000), indicating the model’s lower splashing efficiency, even with ϕ =0. Incidentally, Anderson and Haff’s (1991) numerical experiments for⃒⃒uniform sand grains of 230 μm yielded Nspl ≈ 1.5⃒uimp ⃒ (three times
larger than the results of our Monte Carlo simulations), and Gauer(2001) used this relationship to parameterize the number of splashedparticles. Therefore, we attempted a simulation using three times the
Nspl value obtained for ϕ = 0. The resulting mass transport rate during
the initial 30 s at u*f = 0.3 m s− 1 is shown by an asterisk in Fig. 4a, and isclose to the formula derived by Nishimura and Hunt (2000). However,these results greatly exceed the field observations shown in Fig. 4b.
Hence, further study is needed on the model’s reproducibility for6

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

0.3 m s− 1 are compared with the data observed in Antarctica by Nishi­mura and Nemoto (2005) under a similar wind condition of u*f = 0.28 m
s− 1. Considering that the conditions of the snowpack at the time ofobservation are unknown, and that the measurement uncertainty may
be large, especially at lower heights where the vertical gradient of thedensity and velocity of particles is large, the modelled fluxes are in
reasonable agreement with the observations, both in magnitude andprofile shape. It is notable in Fig. 5a that the reduction in the verticalgradient of flux is reproduced above a height of 0.1–0.2 m, similar to thenumerical simulations of Nemoto and Nishimura (2004), Nemoto et al.
(2004), and Melo et al. (2022). This kink in the profile reflects thetransition from saltation to suspension, with the mass fluxes in thesaltation layer being orders of magnitude larger than those aloft in thesuspension layer. As confirmed by previous field measurements (e.g.,
Takeuchi, 1980; Nishimura and Nemoto, 2005; Nishimura et al., 2014),
the vertical profile of the horizontal mass flux of saltating snow particles
can be approximated by()zq(z) = q0 exp −
,
(17)
Lq

wind-tunnel experiments at low wind speeds.
The mass transport rate varies significantly with time, especially
under high-wind-speed conditions. The transport rate increases rapidly
after the onset of drifting and decreases gradually as the near-surface
flow slows down because of the particle drag, eventually approachingequilibrium. The model results shown by filled symbols in Fig. 4a are theequilibrium transport rates obtained after the flow and particle statistics
reached the steady state throughout the computational domain. There isa significant decrease from the initial transport rates (average of the first30 s), as indicated by the error bars. These overshoot phenomena havebeen predicted numerically by Anderson and Haff (1991) and confirmed
experimentally by Shao and Raupach (1992), who reported that theminimum distance for saltation to reach equilibrium is approximately
15 m. Even longer fetches (>60 m) were required for the snow transport
rate to reach equilibrium in the two-dimensional numerical simulations
reported by Zhang and Zhou (2023). Thus, wind-tunnel experiments
generally measure the transport rates at some point in the early stages ofsaltation growth. The model results are therefore in reasonable corre­
spondence with the wind-tunnel measurements, at least for large u*f .
Fig. 4b presents a comparison with field data observed on flat snow
surfaces with sufficiently long fetches. The data of Kobayashi et al.
(1970) were observed 700 m downwind from the edge of a homoge­
neous snow surface and were generally obtained under low values of thesurface hardness, except for the two smallest values of Q, for which thedata were obtained under hard surface conditions. Takeuchi et al.
(1975) and Takeuchi (1980) presented data observed in the saturation
zone (350 m downwind from the edge), where snowdrift accumulation
had just begun to develop in the leeward direction, and data observed in
the neutral zone, where neither erosion nor deposition were observed.
As shown in the figure, the model results are within the scatter of theobserved data. Under low-wind-speed conditions, the results for the
non-cohesive case are consistent with Kobayashi’s data from the loose
surface and tend to approach the hard surface data as the cohesive en­ergy increases. Under high-wind-speed conditions, the three cases ofdifferent cohesive energies tend to converge on a single curve, which aremost consistent with Takeuchi’s data from the saturation zone, where
the drifting snow and the air flow are likely to be in equilibrium.

where q0 is the reference mass flux and Lq is the characteristic decay
length. The upper height limit for which Eq. (17) fits to the measured
flux has been used as an estimate of the saltation layer height (Sato et al.,
2001; Okaze et al., 2012). The fitted results for the fluxes simulated with
ϕ = 3 × 10− 9 J are shown in Fig. 5b, in which the height hq (determinedas the height where the deviation of Eq. (17) exceeds 10% of the flux) is
indicated by arrows. The values of q0 , Lq , and hq obtained from thedifferent u*f cases are listed in Table 2. The height hq increases with u*f atsmaller velocities, but saturates (or even decreases slightly) at larger u*f .
The length Lq increases monotonically with u*f , indicating that thevertical gradient of the flux (and hence the difference between thesaltation and suspension fluxes) decreases with u*f , consistent with theobservations of Nishimura and Nemoto (2005). The saltation flux esti­
mated as q0 Lq contributes more than 96% of the total mass transport
when u*f = 0.3 m s− 1, but the contribution decreases to 82% in the caseof u*f = 0.8 m s− 1. The relative increase in the contribution of the sus­
pension layer causes an early deviation from the saltation-flux profile(Eq. (17)), which is likely to lead to the saturation of hq .

3.2.2. Vertical profile of horizontal mass flux
Fig. 5 shows vertical profiles of the mean horizontal mass flux q(z). In
Fig. 5a, the fluxes simulated using different cohesive energies at u*f =Fig. 5. Vertical profiles of the mean horizontal mass fluxes: (a) fluxes at u*f = 0.3 m s− 1 modelled with different cohesive energies compared with fluxes observed in
Antarctica by Nishimura and Nemoto (2005) (u*f = 0.28 m s− 1); (b) fluxes at different u*f modelled with ϕ = 3 × 10− 9 J. In (b), dashed lines represent Eq. (17) fitted
to the simulated fluxes within the range 0.02 ≤ z ≤ 0.5 m, the arrows indicate the height hq , determined as the lowest height for which the deviation exceeds 10% of
the flux, and horizontal solid lines indicate the selected heights for which the distribution of particle diameters are shown in Fig. 6.
7

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

Fig. 6. Probability density function (p.d.f.) of airborne particle diameters for different values of u*f : (a) 0.3, (b) 0.4, (c) 0.6, and (d) 0.8 m s− 1. Solid lines represent all
particles, including both aerodynamically entrained and splashed particles, while short-dashed lines represent only particles originating from the aerodynamic
entrainment. The long-dashed line denotes the p.d.f. prescribed for the snow-bed particle diameters.

experimental and observational results (Sugiura et al., 1998; Nishimura
and Nemoto, 2005). The bimodal shape of the distribution at interme­
diate heights (z = 0.042, 0.088 m), especially prominent at lower windspeeds, reflects the contribution of large particles saltating (withrepeated rebounding) on the surface, which is also common to previous
simulation results (Nemoto and Nishimura, 2004; Melo et al., 2022). The
transition probability profiles shown later (Fig. 8) indicate that particles
drifting from various heights remain in a similar height range for a
relatively long time. With respect to the particle feeding process, thecontribution of splash entrainment is generally dominant, except for thelowest u*f case of these simulations using ϕ = 3 × 10− 9 J. Additional
simulations in which rebounding at the surface is supressed by setting
Prbd = 0 result in the disappearance of the secondary peaks at large di­ameters, an increase in the contribution of aerodynamic entrainment,
and a significant decrease in the total mass transport rate (figure notshown). These results confirm the importance of a chain reaction inwhich particles saltating with repeated rebounding over the surface
splash off more particles from the snow bed.

Table 2
Parameters relevant to the saltation layer. q0 and Lq are the parameters of the
exponential profile of the mass flux (Eq. (17)), while hq and hF are the height ofthe saltation layer estimated based on the fit of Eq. (17) and the balance of forces
acting on individual airborne particles, respectively.
u*f [m s− 1]

q0 [kg m− 2 s− 1]

Lq [cm]hq [cm]hF [cm]

0.3
0.35
0.4
0.5
0.6
0.8

6.4 × 10− 2
0.19
0.43
1.1
1.9
3.5

1.6
2.1
2.4
2.8
3.2
4.0

8.3
11.5
12.5
11.5
11.1
11.7

9.1
10.8
12.3
13.9
14.8
16.5

3.2.3. Particle size distribution
Fig. 6 shows the probability density function (p.d.f.) of the airborne
particle diameters at different heights (indicated in Fig. 5b). In Fig. 6,
the solid lines represent all particles, including those both aero­dynamically entrained and splashed from the snow bed, while the shortdashed lines represent aerodynamically entrained particles only. Here­after, the results correspond to simulations using the standard value of ϕ
= 3 × 10− 9 J.
As expected, the particle size distribution is similar to that of thesnow bed at lower heights, the distribution shifts to smaller diameters as
the height increases, and the distribution at higher positions extends tolarger diameters as the wind speed increases. These changes in the p.d.f.
with height and wind speed are qualitatively consistent with previous

3.2.4. Net force acting on drifting particles
In the Lagrangian framework, the acceleration of individual particles
at different heights can be assessed directly by examining the net(aerodynamic and gravitational) force acting on the individual particles.
The aerodynamic and gravitational forces are evaluated by multiplying
the first and second terms on the right-hand side of Eq. (7a) by theparticle mass, respectively. Fig. 7 shows the vertical profiles of the net
force normalized by the gravitational force acting on each particle,
8

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

Fig. 7. Vertical profiles of the mean net force normalized by the gravitational force acting on individual particles, averaged over all particles in logarithmicallyspaced height ranges: (a) profiles over the entire height range; (b) enlarged linear profiles in the dashed rectangular region in (a). Fx and Fz are the streamwise andvertical components of the aerodynamic force, respectively, and Fg is the gravitational force. Data are smoothed using a 1-2-1 vertical filter. Arrows in (b) indicate theheight hF , determined as the lowest height at which the magnitude of the vertical component of the normalized net force decreases below 0.1.

Fig. 8. Profiles of the vertical transition probability of particles originating from different height ranges over the time interval Δttr = 1 s. Individual panels are fordifferent values of u*f : (a) 0.35, (b) 0.5, and (c) 0.8 m s− 1. The numbers associated with each curve denote the total probabilities (%) of particles remaining in the air
after 1 s.

averaged over all particles within logarithmically spaced height ranges.
The figure indicates that the particles at heights below 0.1–0.2 m are, on
average, accelerated in the streamwise direction and decelerated in thevertical direction, thereby producing the streamwise-elongated para­
bolic trajectories observed in the saltation layer (e.g., Sugiura andMaeno, 2000). At higher levels, the gravitational force is almost coun­teracted by the vertical component of the aerodynamic force and thestreamwise aerodynamic force is generally small (though not zero). This
force balance implies that the particles at those heights are almost sus­pended in the air and follow the flow with only marginal acceleration.
Therefore, the lower height limit hF above which the net force becomes
sufficiently small provides another estimate of the saltation layer height.
The arrows in Fig. 7b show hF , estimated as the height at which themagnitude of the normalized net force in the vertical direction decreases
below 0.1 (i.e., 90% of gravity is counteracted by the aerodynamic
force). The values of hF are also listed in Table 2. In contrast to hq , theheight hF increases monotonically with u*f . The streamwise componentof the aerodynamic force near hF increases with wind speed. These facts
imply that, under higher wind speeds, larger particles jump high, but areforced to return to the surface while being accelerated in the streamwise
direction, forming higher, more elongated parabolic trajectories. Thus,
the mean saltation height, estimated purely from the motion of indi­vidual particles, increases with wind speed. It should be noted that anarbitrary threshold of 90% is used here because the net force is nevercompletely zero at any altitude as long as the gravity is present.
Therefore, the value of hF should not be taken as an exact value for the9

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

saltation layer height; however, the analysis in the next subsection
confirms that hF gives a reasonable estimate for the transition heightfrom saltation to suspension.

3.3. Horizontal structure of the drifting snow
3.3.1. Instantaneous correspondence with the near-surface flow
As often experienced (and as shown in Fig. 3 and the animation
provided as Supplementary Video 2), snow drifting over a flat surface
forms streamwise-elongated structures, which “merge, bifurcate, andmeander laterally back and forth as they move downwind” (Bass, 2008).
These are called “snow streamers” or “snow snakes”. The spatial rela­tionship between snow streamers and the instantaneous flow field is
exemplified in Fig. 9, in which the horizontal positions of large number
densities calculated including all particles present within the height
range 0.01 < z ≤ 0.1 m in the saltation layer are marked on the hori­zontal map of instantaneous flow velocities at the lowest grid level (z =0.025 m). The particle number density ρn [m− 3] at a given height follows
a roughly log-normal distribution (figure not shown). Thus, the plotting

3.2.5. Vertical transition probability of drifting particles
The mean vertical displacement over time of particles at each height
can be expressed by the transition probability:
n
∑

(

)

Ptr Δttr , z, zorg =

j=1

(
)
Nptr tj + Δttr , z|tj , zorg
n
∑
j=1

(
)
Np tj , zorg

,

(18)

where Nptr is the number of particles moving from the original heightzorg to height z during the time interval Δttr (=1 s), Np is the total numberof particles present at zorg at time tj when the j-th dataset is saved (seeSection 2.3), and n is the number of datasets. The vertical profiles of Ptr
for different ranges of zorg are shown in Fig. 8, illustrating the differencein particle fates depending on their original positions. After an intervalof 1 s, most particles at levels greater than the estimated height of thesaltation layer (zorg > 0.15 m) remain suspended in the air, and exhibit atendency to disperse further upward. As the wind speed increases, theprobability distribution of these particles extends not only upward, butalso to lower heights, and the total probability decreases. This is because
larger particles, which are more abundant in the air under higher windspeeds (Fig. 6), fall faster, as implied by the greater negative verticalforce balance shown in Fig. 7. Most particles at heights close to thesurface (zorg < 0.07 m) settle to the surface within 1 s, while theremaining particles continue saltating at a height of 0.01–0.04 m, pro­ducing the minor peak of Ptr , with a very small percentage of particles
translated to higher levels. This minor peak of Ptr disappears in therebound-inhibited simulations (figure not shown), indicating that it iscomposed of particles that impacted the surface and rebounded orrepeated the process. At intermediate levels near the saltation layerheight (0.07 < zorg < 0.15 m), approximately 60% of the particles residein the air after 1 s and tend to displace upward, which represents thetransition of particle motion from the saltation mode to the suspension
mode (Bagnold, 1941).

criterion is set as (logρn )’ > σ LN , where (logρn )’ represents a deviation oflogρn from its instantaneous horizontal mean and σLN is the standard

deviation of logρn . Hereafter, a prime (ʹ) indicates deviation from theinstantaneous horizontal mean. The correspondence between thestreamers and the instantaneous flow field is remarkable, with rows
containing multiple marks concentrated in the narrow areas of updrafts(wʹ > 0) with slow streamwise velocities (uʹ < 0). Those areas adjacent

to or enclosed by a larger-scale longitudinal band of fast streamwise
velocities (e.g., 2.5 < y < 5.5 m) are particularly favourable for theformation of streamers. In contrast, in a slower-velocity band (e.g., y <2.5 m), the marks are sparser, despite a comparable number of updraft
regions.
3.3.2. Quadrant analysis

Fig. 10 shows the joint p.d.f. of the turbulent velocity components uʹ
and wʹ, calculated using all datasets for u = 0.5 m s− 1. In this figure, the*fgrey shading with thin dashed line contours represents the joint p.d.f. ofthe velocities obtained from all grid points at the lowest level. As ex­pected, in the near-wall boundary layers (Wallace, 2016), uʹ and wʹ

sampled from the lowest grid points are negatively correlated with eachother, constituting the negative Reynolds stress at that level, and their p.
d.f. is largely distributed in the quadrants of ejection (uʹ < 0, wʹ > 0) andsweep (uʹ > 0, wʹ < 0). The overlaid bold line contours in Fig. 10a
represent a similar p.d.f., but for grid points where the particle number

Fig. 9. Instantaneous distribution of dense particle clouds detected in the saltation layer (square dots) plotted over the coloured distribution of the perturbation flow
velocities at the lowest grid level: (a) streamwise component and (b) vertical component (u*f = 0.5 m s− 1). Dashed lines indicate a large-scale band of faststreamwise velocities.
10T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

Fig. 10. Joint p.d.f. of the turbulent velocity components uʹ and wʹ obtained at the lowest grid level in the case of u*f = 0.5 m s− 1. σ u and σ w are the standard
deviations of uʹ and wʹ, respectively. The grey shading and thin dashed line contours represent the joint p.d.f. of all data at the same level. In (a), the solid line

contours denote the joint p.d.f. calculated only from data at grid points where the particle number density within 0.01 < z ≤ 0.1 m satisfies (log ρn )ʹ > σLN . In (b), thesolid and dashed line contours denote the joint p.d.f. from data at grid points where the number densities within 0 < z ≤ 0.01 m of the freshly splashed particles andfreshly entrained particles exceed the threshold, respectively. The line contours are drawn for 25, 50, and 75% of the maxima.

density in the height range 0.01 < z ≤ 0.1 m is sufficient to satisfy(or streamers), which are associated with ejection motions, are not a
direct result of local enhancement of the particle supply at the same
locations. Rather, they are formed when airborne particles near thesurface, regardless of their origin, gather in narrow regions under theeffects of the flow very close to the surface.

(logρn )’ > σLN . Compared with the p.d.f. for the entire flow, the joint p.
d.f. for velocities at these selected grid points is skewed towards theejection quadrant, indicating that strong ejection events are involved in
the formation of particularly dense particle clouds, which constitute thesnow streamers seen in Fig. 9. The skewness towards the ejection issimilar to that found by Berk and Coletti (2020) from their wind-tunnel
experiments using size-selected glass spheres. To investigate whether
there is any correspondence between the development of dense particle
clouds and the particle feeding process from the snow bed, a similar
analysis is performed for the number densities calculated using only
“freshly entrained” and “freshly splashed” particles, i.e., particles thatwere entrained and splashed at their most recent occurrence (with in­tervals Δtae and Δtspl ), respectively, and have not yet experienced arebound. For this analysis, the respective number densities ρne and ρns
are calculated over a height range 0 < z ≤ 0.01 m, which is below therange used to calculate ρn . The subscripts e and s denote the freshly
entrained and freshly splashed particles, respectively. The results areshown in Fig. 10b, where the bold contours represent the joint p.d.f. of()’
velocities obtained at grid points where logρne,s > σLNe,s is satisfied.
The figure clearly shows that aerodynamic entrainment is associatedwith the most energetic sweep motions. This is because only the strong
stresses caused intermittently by these motions exceed the entrainment
threshold, ρu2*t , in the equilibrium state where the surface stresses are(on average) reduced from the freestream value by the drag of airborne
particles. Consequently, the contribution of aerodynamic entrainment to
the total particle supply is only marginal, as already shown by Fig. 6; thisis also evidenced through this analysis by the mean number density of
the freshly entrained particles (ρne ∼ 105 m− 3) being two orders ofmagnitude smaller than that of the freshly splashed particles (ρns ∼ 107m− 3). Although the p.d.f. of velocities related to the large densities offreshly splashed particles (solid line contours in Fig. 10b) is closer to thebackground p.d.f. of the entire flow, the distribution is still skewed to­wards the sweep quadrant, because the faster flow of the sweep motions
accelerates the particles impacting the surface more efficiently. The
skewness of the p.d.f. for the splash entrainment tends to be smaller for
higher u*f cases, but does not show a preference for the ejection quad­rant in any scenario (figure not shown). Thus, the dense particle clouds

3.3.3. Composite of the flow inducing dense particle clouds near the surface
The mean spatial structure of the flow causing the dense particle
clouds near the surface can be deduced by compositing the conditionallysampled flow variables as)∑n ∑mj ʹ( ii)(i=1 φ xr + rx , yr + ry , z
j=1
∑n
,
(19)
〈φʹ〉 rx , ry , z, zr =j=1 mjwhere 〈φʹ〉 denotes an ensemble average of any flow variable φʹ sampled() ()
i = 1, ⋯, mj , for which thearound detection position xir , yir , zr
detection condition is satisfied in the j-th dataset, rx and ry are thestreamwise and lateral displacements from the detection position,
respectively, n is the total number of datasets, and mj is the number ofdetected points in the j-th dataset. In this analysis, to avoid redundantsampling of the same particle cloud, the detection condition for eachdataset is (logρn )’ > 2σ LN , while ρn is calculated over a height range 0.01
< z ≤ 0.1 m, as above.
Fig. 11 shows the horizontal distributions of the composited flow
field sliced at z = 0.025 m. The figure implies that the dense particle
clouds generally occur near the upstream ends of internal ejection re­gions (uʹ < 0 and wʹ > 0), where near-surface lateral flows converge

towards the low-pressure area. The composited low-speed (uʹ < 0) re­gion is elongated in the streamwise direction, measuring about 0.9 m(length) × 0.1 m (width). This is longer than the average streamer length(0.5 m) estimated by Huang and Wang (2016) from numerical simula­tions in a similar setting to that employed in our study, suggesting that
dense particle clouds are generally encompassed within the flow struc­ture that produces them. The streamwise–vertical section at ry = 0
shown in Fig. 12 indicates that the cores of the dense particle clouds tend
to be located on the leeward side of tall high-pressure zones (Fig. 12c),
which are presumably induced by the stagnation of the streamwise
11

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

Fig. 11. Horizontal cross-sections at the lowest grid level (z = 0.025 m) of the ensemble-averaged flow around the position of particularly dense particle clouds
detected at that level: (a) streamwise velocity component 〈uʹ〉/u*f , (b) cross-streamwise velocity component 〈vʹ〉/u*f , (c) vertical velocity component 〈wʹ〉/ u*f , and (d))(
pressure 〈pʹ〉/ ρu2 , obtained from the case of u = 0.5 m s− 1. Solid and dashed line contours represent positive and negative values, respectively, and are drawn at*f*fintervals of (a) 0.15, (b) 0.1, (c) 0.05, and (d) 0.1.

velocity component near the surface (Fig. 12a). The height of the
low-speed region is confined to within 0.2 m, in contrast to the updraft(wʹ > 0) region, which extends higher up and farther downstream

3.3.4. Relationship with large-scale turbulent motion
The existence of streamwise-elongated low-speed regions near thesurface, known as low-speed streaks, is a ubiquitous feature of near-wallboundary-layer flows (Robinson, 1991). The composited flow systemdescribed in the previous subsection closely resembles the structure of
the near-wall characteristic eddies obtained from direct numericalsimulations of channel flow (Moin and Moser, 1989), such as an ejection
accompanied by a low-speed streak sandwiched by counter-rotating
streamwise vortices that lift downstream, although their spatial scaleis much smaller, as the typical streak spacing (known to be 100 viscouswall units) is 1.5–4 mm under the conditions of this study. The present
analysis indicates that the near-surface streak structure, which is known
to develop through the interaction with small-scale eddies, plays anessential role in the formation of snow streamers, even in high Reynolds
number flows in the atmospheric surface layer.
Another remarkable feature is that the composited low-speed streaks
are enclosed by large high-speed (uʹ > 0) regions that span most of
Figs. 11–13. This reflects the occurrence of streamers being preferen­
tially biased towards the large-scale high-speed bands, as exemplified inFig. 9. Regarding the particle supply, a larger number of particles areseeded to the flow in the large-scale high-speed regions than in their
low-speed counterparts, owing to the enhanced aerodynamic entrain­
ment resulting from larger shear stresses and the increased splashingcaused by the acceleration of impacting particles, as indicated by
Fig. 10b. The entrained and splashed particles are then gathered by thenear-surface convergent flow to form dense particle clouds in the in­ternal low-speed streaks. The activity of near-surface eddies associatedwith the low-speed streaks is modulated by the coherent turbulent(Fig. 12b). These vertical distributions are consistent with thewind-tunnel results of Berk and Coletti (2020), who showed that parti­
cles are concentrated in high strain regions of the near-surface turbulent
flow. The cross-streamwise–vertical sections sliced at different stream­
wise positions shown in Fig. 13 indicate that the internal low-speed
region is situated between a pair of regions of mean streamwise
vorticity with opposite signs. Near the surface, lateral flows from thehigh-speed regions on both sides converge, inducing the updraft in the
middle, and diverge sideways above the vorticity cores. The centre of thevorticity pair lifts slightly downstream. It should be noted, however, that
this ensemble-averaged result does not imply the dominance ofcounter-rotating vortex pairs in the near-surface region. In individual
realizations, small-scale single vortices are more commonly identified
near the boundary between the inner low-speed region and the sur­rounding high-speed regions (figure not shown).
With respect to the formation mechanism of snow streamers, Huang
and Wang (2016) focused on the suction of small-scale vortices, which
may attract saltating snow particles. In contrast, the results of ouranalysis show that the near-surface velocity patterns produced by thecooperative effects of distributed small-scale vortices, rather than indi­vidual vortices themselves, are more closely correlated with the snow
streamers. This view does not contradict the finding of Huang and Wang
(2016) that snow streamers form in regions flanked on both sides bylarge-magnitude vorticity (see Fig. 8b in their paper).

12

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

Fig. 12. Streamwise–vertical cross-sections at ry = 0 of the ensemble-averaged flow around the position of particularly dense particle clouds detected at the lowest)(
grid level: (a) streamwise velocity component 〈uʹ〉/u , (b) vertical velocity component 〈wʹ〉/u , and (c) pressure 〈pʹ〉/ ρu2 , obtained from the case of u = 0.5 m
*f

*f

*f

*f

s− 1. Solid and dashed line contours represent positive and negative values, respectively, and are drawn at intervals of (a) 0.2, (b) 0.1, and (c) 0.2. The grey line
contours in (b) and (c) replicate the contours of 〈uʹ〉/u shown in (a).
*fmotion of much larger spatial scales—the small-scale velocity fluctua­
tions near the surface are amplified within the high-speed regions of thelarger coherent motions (Marusic et al., 2010). As demonstrated byAndreolli et al. (2023), the nonlinear energy transfer between large- andsmall-scale motions, most likely through the stretching of small eddies,is involved in the modulation mechanism. Besides the direct conse­
quence of the more intense internal near-surface flow system, theenhanced activity of small-scale eddies causes the meandering motionsof the internal low-speed streaks. Both effects lead to the intensificationand fluctuation of the snow streamers formed in the large-scale high-­speed regions.
In the atmospheric boundary layers, very large eddies impinging
from aloft generate wide high-speed regions (or very large-scale sweeps)
near the surface, causing very high spatio-temporal intermittency in thenear-surface velocity fluctuations (Hunt and Morrison, 2000). Sucheddies are generally scaled with the height of the outer or internalboundary layers and/or regulated by nearby topography. Field obser­
vations have shown that the velocity fluctuations associated with suchlarge eddies have a significant impact on the concentration and flux ofdrifting particles. Aksamit and Pomeroy (2018a,b) observed that theinitiation of intensive blowing-snow transport and the increase of theconcentration and particle number flux near the surface corresponded tothe passage of low-frequency high-speed motions. They also highlighted
the role of these high-speed motions that amplify the high-frequency
fluctuations in blowing snow density by modulating high-frequency
velocity fluctuations. Bass and Sherman (2005) identified that
extremely large transport of saltating sand particles occurs in the form of
ensembled streamers (similar to that seen in the high-speed band in
Fig. 9) during the periods of intermittent high-speed airflow. The latter
authors estimated that the spatial scale of sand streamers is determinedby the eddies scaled with the height of the internal boundary layer
developing over the observation site. Although our simulations do notinclude the influence of such boundary-layer scale eddies, the results are
consistent with these observations in the sense that the large-scale
high-speed motions increase the particle supply from the surface andactivate the small-scale eddies that contribute to the formation of indi­vidual streamers.
4. Conclusions
We have described a Lagrangian particle dispersion model coupledwith large-eddy simulation code based on the central-moment lattice
Boltzmann method. This model was used to simulate drifting snow in aneutrally stratified turbulent flow within the atmospheric surface layer.
The model reproduces typical features of drifting snow observed in thefield, such as the dependency of the total mass transport rate on the flow
velocity, the kink in the vertical profile of the mass flux near the salta­tion layer height, and the variations in particle size distribution withflow velocity and height from the snow bed.
The saltation layer height, determined directly from the inspection offorces acting on the airborne particles, increases monotonically with theflow velocity (or the shear stress) above the saltation layer. This isdifferent from conventional estimates, which tend to saturate withincreasing flow velocity. The vertical transition probability, which wasintroduced to evaluate the mean tendency over time of the verticaldisplacement of individual particles at each height, supports the gradualtransition from saltation to suspension around the saltation layer height
estimated by the present method.
Regarding the horizontal structure of the saltation layer, snow
streamers were found to be closely associated with the appearance ofsmall-scale low-speed streaks formed in the instantaneous flow fieldnear the surface. Particularly dense particle clouds are most often foundnear the upstream ends of these low-speed streaks, the result of particles
being concentrated by near-surface convergent flows. From a larger13

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

Fig. 13. Cross-streamwise–vertical cross-sections of the ensemble-averaged flow around the position of particularly dense particle clouds detected at the lowest gridlevel sliced at different rx : (a) 0 m, (b) +0.1 m, (c) +0.2 m, (d) +0.3 m, and (e) +0.4 m, obtained from the case of u*f = 0.5 m s− 1. The colour shading indicates 〈uʹ〉/)(u*f , the vectors represent 〈vʹ〉 /u*f , 〈wʹ〉 /u*f , and the solid and dashed line contours drawn at intervals of 0.05 represent positive and negative values of the()
streamwise vorticity 〈ωʹ 〉/ u /Δx , respectively.
x*fscale perspective, the low-speed steaks appear almost everywhere near
the surface, whereas dense particle clouds are selectively associatedwith those streaks that are enclosed by high-speed regions of muchlarger spatial scale. The reasons for this preference are the greater par­ticle supply through the aerodynamic and splash entrainment processes(resulting from the higher flow velocities) and the enhanced activities ofsmall-scale turbulent eddies (most likely owing to increased stretching)in the large-scale high-speed regions. Although the limited computa­
tional domain used in this study prevents generalization, we speculatethat the preferential occurrence of snow streamers, associated with thelarge-eddy-modulated low-speed streaks in the vicinity of the snow
surface, may be one of the fundamental mechanisms behind the spatiotemporal intermittency of the drifting snow in the surface layers. Further
analyses on whether saltating particles affect the near-surface flows and
14

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

their modulation in response to the large-scale coherent motions are left
for future research.
In terms of model developments, this study only dealt with an ideal
situation, i.e. drifting of spherical particles over a flat and homogeneous
snow surface under a neutrally stratified atmosphere. Therefore, areasonable extension of the model enabling its application to more
realistic situations would be to include the effects of the buoyancy andsublimation of particles with the aid of appropriate solvers for temper­
ature and humidity in the flow. Consideration of non-spheric particles isalso an issue for future research. In such studies, a model of the dragcoefficient (e.g., Tagliavini et al., 2021) and general parametrizations of
aerodynamic lift forces for particles with complex geometry should be
explored. Moreover, using the adaptability of the lattice Boltzmann
method to complex boundaries, it is feasible to extend the model tosimulate the time-varying surface topographies created by erosion anddeposition (around obstacles or on flat surfaces), and the resulting
feedback to the turbulent flow. Possible applications of the model
framework are not limited to the dispersion of particulate matter such as
dust, pollen, sprays, and other aerosols, but also include Lagrangian
trajectories of any passive scalar quantities represented as masslessparticles.

editing, Visualization, Software, Methodology, Formal analysis.
Masayuki Kawashima: Writing – review & editing, Methodology. Kou
Shimoyama: Writing – review & editing, Methodology. Naoyuki
Onodera: Writing – review & editing, Software, Methodology. Yuta
Hasegawa: Writing – review & editing, Methodology. Atsushi Inagaki:
Writing – review & editing, Project administration, Methodology.
Declaration of competing interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.
Data availability
Data will be made available on request.
Acknowledgements
The authors thank Hirofumi Niiya and Takane Matsumoto (Niigata
University) for their fruitful discussions, and Kazuya Ono (Technical
Division, Institute of Low Temperature Science, Hokkaido University)
for his support with the server maintenance. This work was supported inpart by JSPS KAKENHI (Grant Number JP20K04057) and by the Joint
Research Program of the Institute of Low Temperature Science, Hok­kaido University (Grant Numbers 22G004, 23G008). The authors also
thank Edanz (jp.edanz.com) for English language editing.

CRediT authorship contribution statement
Tsutomu Watanabe: Writing – original draft, Visualization, Super­
vision, Software, Resources, Methodology, Funding acquisition, Formal
analysis, Conceptualization. Shuhei Ishikawa: Writing – review &Appendix A. Subgrid-scale velocity experienced by particles
In the calculation of the particle trajectory, not only resolved fluid velocities but also unresolved, subgrid-scale (SGS) velocities are considered (Eq.
(7)). The SGS velocities experienced by the particles are estimated based on the Lagrangian stochastic model proposed by Weil et al. (2004) with somemodifications as√̅̅̅̅̅̅̅̅)(− us
1 1 dσ 2s
2σ2s
2
dus =
dξ,
(20)
dt +
us + ∇σs dt +2 σ2s dtτs
τswhere τs is the timescale of the SGS velocity, σ 2s is the SGS velocity variance evaluated in the large-eddy simulation model, and dξ is a vector composedof Gaussian white noises ξi ∈ N(0, dt). Following Knorps and Pozorski (2021), τs is parametrized asτs = aτmin(Δx, z + z0s )
,
2σ s

(21)

where aτ , which is a parameter representing the effect of the particle inertia, is assumed constant (=0.5) considering the large Stokes numbers in thenear-wall region of our simulations (Table 1). For the particles below the lowest grid level zb , σs is modified as follows:

σ s (z) = σ s (zb )f1 (z) + σ nw f2 (z) (0 < z < zb ),(22)

where σnw is the near-wall standard deviation of velocities estimated as σnw = 1.7u based on DNS results of channel flow (e.g., Moser et al., 1999) andf1 and f2 are logarithmic interpolation functions defined as
f1 (z) =ln(z + z0s ) − ln(z0s ),ln(zb + z0s ) − ln(z0s )(23a)

f2 (z) =ln(zb + z0s ) − ln(z + z0s ).
ln(zb + z0s ) − ln(z0s )(23b)

For the levels very close to the surface, z+(zu/ν)<20 , a linear dumping of z+ /20 is applied to the value calculated from Eq. (22).
Appendix B. Supplementary data
Supplementary data to this article can be found online at https://doi.org/10.1016/j.jweia.2024.105783.

15

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783

T. Watanabe et al.

References

Gromke, C., Manes, C., Walter, B., Lehning, M., Guala, M., 2011. Aerodynamic roughness
length of fresh snow. Boundary-Layer Meteorol. 141, 21–34. https://doi.org/
10.1007/s10546-011-9623-3.
Groot Zwaaftink, C.D., Diebold, M., Horender, S., Overney, J., Lieberherr, G.,
Parlange, M.B., Lehning, M., 2014. Modelling small-scale drifting snow with a
Lagrangian stochastic model based on large-eddy simulations. Boundary-Layer
Meteorol. 153, 117–139. https://doi.org/10.1007/s10546-014-9934-2.
Han, M., Ooka, R., Kikumoto, H., 2020. Validation of lattice Boltzmann method-basedlarge-eddy simulation applied to wind flow around single 1:1:2 building model.
J. Wind Eng. Ind. Aerod. 206, 104277 https://doi.org/10.1016/j.
jweia.2020.104277.
Huang, N., Wang, Z.S., 2016. The formation of snow streamers in the turbulent
atmosphere boundary layer. Aeolian Res 23, 1–10. https://doi.org/10.1016/j.
aeolia.2016.09.002.
Hunt, J.C.R., Morrison, J.F., 2000. Eddy structure in turbulent boundary layers. Eur. J.
Mech. B 19, 673–694. https://doi.org/10.1016/S0997-7546(00)00129-1.
Inagaki, A., Kanda, M., Ahmad, N.H., Yagi, A., Onodera, N., Aoki, T., 2017. A numerical
study of turbulence statistics and the structure of a spatially-developing boundary
layer over a realistic urban geometry. Boundary-Layer Meteorol. 164, 161–181.
https://doi.org/10.1007/s10546-017-0249-y.
Islam, A., Gaylard, A., Thornber, B., 2017. A detailed statistical study of unsteady wake
dynamics from automotive bluff bodies. J. Wind Eng. Ind. Aerod. 171, 161–177.
https://doi.org/10.1016/j.jweia.2017.09.009.
Kareem, A., 2020. Emerging frontiers in wind engineering: computing, stochastics,
machine learning and beyond. J. Wind Eng. Ind. Aerod. 206, 104320 https://doi.
org/10.1016/j.jweia.2020.104320.
Knorps, M., Pozorski, J., 2021. Stochastic modeling for subgrid-scale particle dispersion
in large-eddy simulation of inhomogeneous turbulence. Phys. Fluids 33, 043323.
https://doi.org/10.1063/5.0046320.
Kobayashi, D., Kobayashi, S., Ishikawa, N., 1970. Measurement of snow drift using
parallel trenches. Low Temp. Sci. 27, 99–106 (in Japanese with English summary).
http://hdl.handle.net/2115/18104.
Kobayashi, H., 2005. The subgrid-scale models based on coherent structures for rotatinghomogeneous turbulence and turbulent channel flow. Phys. Fluids 17, 045104.
https://doi.org/10.1063/1.1874212.
Kok, J.F., Renno, N.O., 2009. A comprehensive numerical model of steady state saltation(COMSALT). J. Geophys. Res. 114, D17204 https://doi.org/10.1029/
2009JD011702.
Kondo, J., Yamazawa, H., 1986. Bulk transfer coefficient over a snow surface. BoundaryLayer Meteorol. 34, 123–135. https://doi.org/10.1007/BF00120912.
Kosugi, K., Nishimura, K., Maeno, N., 1995. Studies on the dynamics of saltation in
drifting snow. Rep. Natl. Res. Inst. Earth Sci. Disaster Prev. 54, 111–154. https://doi.
org/10.24732/nied.00001063.
Krüger, T., Kusumaatmaja, H., Kuzmin, A., Shardt, O., Silva, G., Viggen, E.M., 2017. The
Lattice Boltzmann Method: Principles and Practice. Springer.
Lenz, S., Schönherr, M., Geier, M., Krafczyk, M., Pasquali, A., Christen, A., Giometto, M.,
2019. Towards real-time simulation of turbulent air flow over a resolved urban
canopy using the cumulant lattice Boltzmann method on a GPGPU. J. Wind Eng. Ind.
Aerod. 189, 151–162. https://doi.org/10.1016/j.jweia.2019.03.012.
Liston, G.E., Sturm, M., 1998. A snow-transport model for complex terrain. J. Glaciol. 44,
498–516. https://doi.org/10.3189/S0022143000002021.
Liston, G.E., Haehnel, R.B., Strum, M., Hiemstra, C.A., Berezovskaya, S., Tabler, R.D.,
2007. Simulating complex snow distributions in windy environments using
SnowTran-3D. J. Glaciol. 53, 241–256. https://doi.org/10.3189/
172756507782202865.
Liu, Z., Zhou, X., Gu, M., 2022. A numerical model for snow drifting simulations on flat
roofs using Lagrangian approach. J. Wind Eng. Ind. Aerod. 222, 104922 https://doi.
org/10.1016/j.jweia.2022.104922.
Marusic, I., Mathis, R., Hutchins, N., 2010. Predictive model for wall-bounded turbulent
flow. Science 329, 193–196. https://doi.org/10.1126/science.1188765.
Melo, D.B., Sharma, V., Comola, F., Sigmund, A., Lehning, M., 2022. Modeling snow
saltation: the effect of grain size and interparticle cohesion. J. Geophys. Res. Atmos.
127, e2021JD035260 https://doi.org/10.1029/2021JD035260.
Moin, P., Moser, R.D., 1989. Characteristic-eddy decomposition of turbulence in achannel. J. Fluid Mech. 200, 471–509. https://doi.org/10.1017/
S0022112089000741.
Moser, R.D., Kim, J., Mansour, N.N., 1999. Direct numerical simulation of turbulent
channel flow up to Reτ = 590. Phys. Fluids 11, 943–945. https://doi.org/10.1063/
1.869966.
Mott, R., Vionnet, V., Grünewald, T., 2018. The seasonal snow cover dynamics: review
on wind-driven coupling processes. Front. Earth Sci. 6, 197. https://doi.org/
10.3389/feart.2018.00197.
Nemoto, M., Nishimura, K., 2004. Numerical simulation of snow saltation andsuspension in a turbulent boundary layer. J. Geophys. Res. 109, D18206 https://doi.
org/10.1029/2004JD004657.
Nemoto, M., Nishimura, K., Kobayashi, S., Izumi, K., 2004. Numerical study of the timedevelopment of drifting snow and its relation to the spatial development. Ann.
Glaciol. 38, 343–350. https://doi.org/10.3189/172756404781815202.
Niiya, H., Nishimura, K., 2017. Spatiotemporal structure of aeolian particle transport onflat surface. J. Phys. Soc. Japan 86, 054402. https://doi.org/10.7566/
JPSJ.86.054402.
Niiya, H., Nishimura, K., 2022. Hysteresis and surface shear stresses during snow-particle
aeolian transportation. Boundary-Layer Meteorol. 183, 447–467. https://doi.org/
10.1007/s10546-022-00688-8.

Ahmad, N.H., Inagaki, A., Kanda, M., Onodera, N., Aoki, T., 2017. Large-eddy simulationof the gust index in an urban area using the lattice Boltzmann method. BoundaryLayer Meteorol. 163, 447–467. https://doi.org/10.1007/s10546-017-0233-6.
Aidun, C.K., Clausen, J.R., 2010. Lattice-Boltzmann method for complex flows. Annu.
Rev. Fluid Mech. 42, 439–472. https://doi.org/10.1146/annurev-fluid-121108145519.
Aksamit, N.O., Pomeroy, J.W., 2018a. The effect of coherent structures in theatmospheric surface layer on blowing-snow transport. Boundary-Layer Meteorol.
167, 211–233. https://doi.org/10.1007/s10546-017-0318-2.
Aksamit, N.O., Pomeroy, J.W., 2018b. Scale interactions in turbulence for mountainblowing snow. J. Hydrometeorol. 19, 305–320. https://doi.org/10.1175/JHM-D-170179.1.
Ammi, M., Oger, L., Beladjine, D., Valance, A., 2009. Three-dimensional analysis of thecollision process of a bead on a granular packing. Phys. Rev. E 79, 021305. https://doi.org/10.1103/PhysRevE.79.021305.
Anderson, R.S., Haff, P.K., 1991. Wind modification and bed response during saltation ofsand in air. In: Barndorff-Nielsen, O.E., Willetts, B.B. (Eds.), Aeolian Grain Transport
1. Acta Mech. Supp. 1. Springer, Vienna, pp. 21–51. https://doi.org/10.1007/978-37091-6706-9_2.
Andre, M., Mier-Torrecilla, M., Wüchner, R., 2015. Numerical simulation of wind loads
on a parabolic trough solar collector using lattice Boltzmann and finite element
methods. J. Wind Eng. Ind. Aerod. 146, 185–194. https://doi.org/10.1016/j.
jweia.2015.08.010.
Andreolli, A., Gatti, D., Vinuesa, R., Örlü, R., Schlatter, P., 2023. Separating large-scale
superposition and modulation in turbulent channels. J. Fluid Mech. 958, A37.
https://doi.org/10.1017/jfm.2023.103.
Araoka, K., Maeno, N., 1981. Dynamical behaviors of snow particles in the saltation
layer. Mem. Natl. Inst. Polar Res. - Special Issue 19, 253–263. http://id.nii.ac.jp/1291/00001187.
Bagnold, R.A., 1941. The Physics of Blown Sand and Desert Dunes. Methuen, London.
Bass, A.C.W., 2008. Challenges in aeolian geomorphology: investigating aeolian
streamers. Geomorphology 93, 3–16. https://doi.org/10.1016/j.
geomorph.2006.12.015.
Bass, A.C.W., Sherman, D.J., 2005. Formation and behavior of aeolian streamers.
J. Geophys. Res. 110, F03011 https://doi.org/10.1029/2004JF000270.
Berk, T., Coletti, F., 2020. Transport of inertial particles in high-Reynolds-number
turbulent boundary layers. J. Fluid Mech. 903, A18 https://doi.org/10.1017/
jfm.2020.597.
Beyers, J.H.M., Sundsbø, P.A., Harms, T.M., 2004. Numerical simulation of threedimensional, transient snow drifting around a cube. J. Wind Eng. Ind. Aerod. 92,
725–747. https://doi.org/10.1016/j.jweia.2004.03.011.
Buffa, E., Jacob, J., Sagaut, P., 2021. Lattice-Boltzmann-based large-eddy simulation ofhigh-rise building aerodynamics with inlet turbulence reconstruction. J. Wind Eng.
Ind. Aerod. 212, 104560 https://doi.org/10.1016/j.jweia.2021.104560.
Chen, S., Doolen, G.D., 1998. Lattice Boltzmann method for fluid flows. Annu. Rev. Fluid
Mech. 30, 329–364. https://doi.org/10.1146/annurev.fluid.30.1.329.
Chen, X., Yu, Z., Liu, Z., Cao, R., Zhu, F., 2021. A prediction model for roof snow
distribution during long-period snowfall using Eular-Lagrangian method. Cold Reg.
Sci. Technol. 190, 103347 https://doi.org/10.1016/j.coldregions.2021.103347.
Clifton, A., Lehning, M., 2008. Improvement and validation of a snow saltation modelusing wind tunnel measurements. Earth Surf. Process. Landforms 33, 2156–2173.
https://doi.org/10.1002/esp.1673.
Clifton, A., Rüedi, J.D., Lehning, M., 2006. Snow saltation threshold measurements in adrifting-snow wind tunnel. J. Glaciol. 52, 585–596. https://doi.org/10.3189/
172756506781828430.
Comola, F., Kok, J.F., Gaume, J., Paterna, E., Lehning, M., 2017. Fragmentation of windblown snow crystals. Geophys. Res. Lett. 44, 4195–4203. https://doi.org/10.1002/
2017GL073039.
Comola, F., Lehning, M., 2017. Energy- and momentum-conserving model of splash
entrainment in sand and snow saltation. Geophys. Res. Lett. 44, 1601–1609. https://doi.org/10.1002/2016GL071822.
Doorschot, J.J.J., Lehning, M., Vrouwe, A., 2004. Field measurements of snow-drift
threshold and mass fluxes, and related model simulations. Boundary-Layer Meteorol.
113, 347–368. https://doi.org/10.1007/s10546-004-8659-z.
Dupont, S., Bergametti, G., Marticorena, B., Simoëns, S., 2013. Modeling saltationintermittency. J. Gehophys. Res. Atmos. 118, 7109–7128. https://doi.org/10.1002/
jgrd.50528.
Feng, Y., Miranda-Fuentes, J., Guo, S., Jacob, J., Sagaut, P., 2021. ProLB: a lattice
Boltzmann solver of large-eddy simulation for atmospheric boundary layer flows.
J. Adv. Model. Earth Syst. 13, e2020MS002107 https://doi.org/10.1029/
2020MS002107.
Gauer, P., 2001. Numerical modeling of blowing and drifting snow in Alpine terrain.
J. Glaciol. 47, 97–110. https://doi.org/10.3189/172756501781832476.
Geier, M., Greiner, A., Korvink, J.G., 2006. Cascaded digital lattice Boltzmann automata
for high Reynolds number flow. Phys. Rev. E 73, 066705. https://doi.org/10.1103/
PhysRevE.73.066705.
Geier, M., Schönherr, M., Pasquali, A., Krafczyk, M., 2015. The cumulant lattice
Boltzmann equation in three dimensions: theory and validation. Comput. Math.
Appl. 70, 507–547. https://doi.org/10.1016/j.camwa.2015.05.001.
Gordon, M., Taylor, P.A., 2009. Measurements of blowing snow, Part I: particle shape,size distribution, velocity, and number flux at Churchill, Manitoba, Canada. ColdReg. Sci. Technol. 55, 63–74. https://doi.org/10.1016/j.coldregions.2008.05.001.

16

T. Watanabe et al.

Journal of Wind Engineering & Industrial Aerodynamics 250 (2024) 105783
Takeuchi, M., Ishimoto, K., Nohara, T., 1975. A study of drift snow transport. J. Jpn. Soc.
Snow Ice 37, 114–121. https://doi.org/10.5331/seppyo.37.114 (in Japanese withEnglish summary).
Takeuchi, M., 1980. Vertical profile and horizontal increase of drift-snow transport.
J. Glaciol. 26, 481–492. https://doi.org/10.3189/S0022143000010996.
Tanji, S., Inatsu, M., Okaze, T., 2021. Development of a snowdrift model with the lattice
Boltzmann method. Prog. Earth Planet. Sci. 8, 57. https://doi.org/10.1186/s40645021-00449-0.
Tominaga, Y., Mochida, A., 1999. CFD prediction of flowfield and snowdrift around abuilding complex in a snowy region. J. Wind Eng. Ind. Aerod. 81, 273–282. https://doi.org/10.1016/S0167-6105(99)00023-9.
Tominaga, Y., Mochida, A., Okaze, T., Sato, T., Nemoto, M., Motoyoshi, H., Nakai, S.,
Tsutsumi, T., Otsuki, M., Uematsu, T., Yoshino, H., 2011. Development of a systemfor predicting snow distribution in built-up environments: combining a mesoscalemeteorological model and a CFD model. J. Wind Eng. Ind. Aerod. 99, 460–468.
https://doi.org/10.1016/j.jweia.2010.12.004.
Tominaga, Y., 2018. Computational fluid dynamics simulation of snowdrift aroundbuildings: past achievements and future perspectives. Cold Reg. Sci. Technol. 150,
2–14. https://doi.org/10.1016/j.coldregions.2017.05.004.
Uematsu, T., Nakata, T., Takeuchi, K., Arisawa, Y., Kaneda, Y., 1991. Three-dimensional
numerical simulation of snowdrift. Cold Reg. Sci. Technol. 20, 65–73. https://doi.
org/10.1016/0165-232X(91)90057-N.
Vionnet, V., Martin, E., Masson, V., Guyomarc’h, G., Naaim-Bouvet, F., Prokop, A.,
Durand, Y., Lac, C., 2014. Simulation of wind-induced snow transport andsublimation in alpine terrain using a fully coupled snowpack/atmosphere model.
Cryosphere 8, 395–415. https://doi.org/10.5194/tc-8-395-2014.
Wallace, J.M., 2016. Quadrant analysis in turbulence research: history and evolution.
Annu. Rev. Fluid Mech. 48, 131–158. https://doi.org/10.1146/annurev-fluid122414-034550.
Wang, Z., Huang, N., Pähtz, T., 2019. The effect of turbulence on drifting snow
sublimation. Geophys. Res. Lett. 46, 11568–11575. https://doi.org/10.1029/
2019GL083636.
Wang, Z., Jia, S., 2018. A simulation of a large-scale drifting snowstorm in the turbulent
boundary layer. Cryosphere 12, 3841–3851. https://doi.org/10.5194/tc-12-38412018.
Watanabe, T., 2004. Large-eddy simulation of coherent turbulence structures associatedwith scalar ramps over plant canopies. Boundary-Layer Meteorol. 112, 307–341.
https://doi.org/10.1023/B:BOUN.0000027912.84492.54.
Watanabe, T., 2009. LES study on the structure of coherent eddies inducing predominant
perturbations in velocities in the roughness sublayer over plant canopies.
J. Meteorol. Soc. Japan 87, 39–56. https://doi.org/10.2151/jmsj.87.39.
Watanabe, T., Shimoyama, K., Kawashima, M., Mizoguchi, Y., Inagaki, A., 2020. Largeeddy simulation of neutrally-stratified turbulent flow within and above plant canopy
using the central-moments-based lattice Boltzmann method. Boundary-Layer
Meteorol. 176, 35–60. https://doi.org/10.1007/s10546-020-00519-8.
Watanabe, T., Takagi, M., Shimoyama, K., Kawashima, M., Onodera, N., Inagaki, A.,
2021. Coherent eddies transporting passive scalars through the plant canopy
revealed by large-eddy simulations using the lattice Boltzmann method. BoundaryLayer Meteorol. 181, 39–71. https://doi.org/10.1007/s10546-021-00633-1.
Weil, J.C., Sullivan, P.P., Moeng, C.H., 2004. The use of large-eddy simulations in
Lagrangian particle dispersion models. J. Atmos. Sci. 61, 2877–2887. https://doi.
org/10.1175/JAS-3302.1.
White, F.M., 1974. Viscous Fluid Flow. McGraw-Hill, New York.
Xing, M., He, C., 2013. 3D ejection behavior of different sized particles in the grain-bed
collision process. Geomorphology 187, 94–100. https://doi.org/10.1016/j.
geomorph.2013.01.002.
Zhang, T., Zhou, X., 2023. Numerical simulations of snow saltation over flat terrain:development of snow transport. J. Wind Eng. Ind. Aerod. 204, 105472 https://doi.
org/10.1016/j.jweia.2023.105472.
Zhou, X., Zhang, T., 2023. A review of computational fluid dynamics simulations ofwind-induced snow drifting around obstacles. J. Wind Eng. Ind. Aerod. 234, 105350
https://doi.org/10.1016/j.jweia.2023.105350.
Zhou, X., Zhang, T., Liu, Z., Gu, M., 2023. An Eulerian-Lagrangian snow drifting modelfor gable roofs considering the effect of roof slope. J. Wind Eng. Ind. Aerod. 234,
105334 https://doi.org/10.1016/j.jweia.2023.105334.

Nishimura, K., Yokoyama, C., Ito, Y., Nemoto, M., Naaim-Bouvet, F., Bellot, H., Fujita, K.,
2014. Snow particle speeds in drifting snow. J. Geophys. Res. Atmos. 119,9901–9913. https://doi.org/10.1002/2014JD021686.
Nishimura, K., Hunt, J.C.R., 2000. Saltation and incipient suspension above a flat particle
bed below a turbulent boundary layer. J. Fluid Mech. 417, 77–102. https://doi.org/
10.1017/S0022112000001014.
Nishimura, K., Nemoto, M., 2005. Blowing snow at Mizuho station, Antarctica. Phil.
Trans. R. Soc. A 363, 1647–1662. https://doi.org/10.1098/rsta.2005.1599.
Okaze, T., Mochida, A., Tominaga, Y., Nemoto, M., Sato, T., Sasaki, Y., Ichinohe, K.,
2012. Wind tunnel investigation of drifting snow development in a boundary layer.
J. Wind Eng. Ind. Aerod. 104–106, 532–539. https://doi.org/10.1016/j.
jweia.2012.04.002.
Okaze, T., Niiya, H., Nishimura, K., 2018. Development of a large-eddy simulation
coupled with Lagrangian snow transport model. J. Wind Eng. Ind. Aerod. 183,
35–43. https://doi.org/10.1016/j.jweia.2018.09.027.
Onodera, N., Idomura, Y., Hasegawa, Y., Nakayama, H., Shimokawabe, T., Aoki, T.,
2021. Real-time tracer dispersion simulations in Oklahoma City using the locallymesh-refined lattice Boltzmann method. Boundary-Layer Meteorol. 179, 187–208.
https://doi.org/10.1007/s10546-020-00594-x.
Pomeroy, J.W., Gray, D.M., Landine, P.G., 1993. The Prairie blowing snow model:
characteristics, validation. operation. J. Hydrol. 144, 165–192. https://doi.org/
10.1016/0022-1694(93)90171-5.
Rice, M.A., Willetts, B.B., McEwan, I.K., 1995. An experimental study of multiple grainsize ejecta produced by collisions of saltating grains with a flat bed. Sedimentology
42, 695–706. https://doi.org/10.1111/j.1365-3091.1995.tb00401.x.
Robinson, S.K., 1991. Coherent motions in the turbulent boundary layer. Annu. Rev.
Fluid Mech. 23, 601–639. https://doi.org/10.1146/annurev.fl.23.010191.003125.
Sato, T., Kosugi, K., Mochizuki, S., Nemoto, M., 2008. Wind speed dependences offracture and accumulation of snowflakes on snow surface. Cold Reg. Sci. Technol.
51, 229–239. https://doi.org/10.1016/j.coldregions.2007.05.004.
Sato, T., Kosugi, K., Sato, A., 2001. Saltation-layer structure of drifting snow observed inwind tunnel. Ann. Glaciol. 32, 203–208. https://doi.org/10.3189/
172756401781819184.
Shao, Y., Li, A., 1999. Numerical modelling of saltation in the atmospheric surface layer.
Boundary-Layer Meteorol. 91, 199–225. https://doi.org/10.1023/A:
1001816013475.
Sato, T., Uematsu, T., Nakata, T., Kaneda, Y., 1993. Three dimensional numericalsimulation of snowdrift. J. Wind Eng. Ind. Aerod. 46/47, 741–746. https://doi.org/
10.1016/B978-0-444-81688-7.50082-6.
Shao, Y., Raupach, M.R., 1992. The overshoot and equilibration of saltation. J. Geophys.
Res. 97, 20559–20564. https://doi.org/10.1029/92JD02011.
Sharma, V., Comola, F., Lehning, M., 2018. On the suitability of the Thorpe–Masonmodel for calculating sublimation of saltating snow. Cryosphere 12, 3499–3509.
https://doi.org/10.5194/tc-12-3499-2018.
Sigmund, A., Dujardin, J., Comola, F., Sharma, V., Huwald, H., Melo, D.B., Hirasawa, N.,
Nishimura, K., Lehning, M., 2022. Evidence of strong flux underestimation by bulk
parametrizations during drifting and blowing snow. Boundary-Layer Meteorol. 182,
119–146. https://doi.org/10.1007/s10546-021-00653-x.
Sørensen, M., 2004. On the rate of aeolian sand transport. Geomorphology 59, 53–62.
https://doi.org/10.1016/j.geomorph.2003.09.005.
Sugiura, K., Nishimura, K., Maeno, N., Kimura, T., 1998. Measurements of snow massflux and transport rate at different particle diameters in drifting snow. Cold Reg. Sci.
Technol. 27, 83–89. https://doi.org/10.1016/S0165-232X(98)00002-0.
Sugiura, K., Maeno, N., 2000. Wind-tunnel measurements of restitution coefficients andejection number of snow particles in drifting snow: determination of splash
functions. Boundary-Layer Meteorol. 95, 123–143. https://doi.org/10.1023/A:
1002681026929.
Syms, G.F., 2008. Simulation of simplified-frigate airwakes using a lattice-Boltzmann
method. J. Wind Eng. Ind. Aerod. 96, 1197–1206. https://doi.org/10.1016/j.
jweia.2007.06.040.
Tagliavini, G., McCorquodale, M., Westbrook, C., Corso, P., Krol, Q., Holzner, M., 2021.
Drag coefficient prediction of complex-shaped snow particles falling in air beyond
the Stokes regime. Int. J. Multiphas. Flow 140, 103652. https://doi.org/10.1016/j.
ijmultiphaseflow.2021.103652.

17