Progress in Earth andPlanetary Science(2021) 8:57
Tanji et al. Prog Earth Planet Sci
https://doi.org/10.1186/s40645-021-00449-0

Open Access

RESEARCH ARTICLE

Development of a snowdrift modelwith the lattice Boltzmann method
Seika Tanji1* , Masaru Inatsu2,3 and Tsubasa Okaze4

Abstract
We developed a snowdrift model to evaluate the snowdrift height around snow fences, which are often installed
along roads in snowy, windy locations. The model consisted of the conventional computational fluid dynamics solverthat used the lattice Boltzmann method and a module for calculating the snow particles’ motion and accumulation.
The calculation domain was a half channel with a flat free-slip boundary on the top and a non-slip boundary on thebottom, and an inflow with artificially generated turbulence from one side to the outlet side was imposed. In addition
to the reference experiment with no fence, experiments were set up with a two-dimensional and a three-dimensional
fence normal to the dominant wind direction in the channel center. The estimated wind flow over the two-dimensional fence was characterized by a swirling eddy in the cross section, whereas the wind flow in the three-dimensional fence experiment was horizontally diffluent with a dipole vortex pair on the leeward side of the fence. Almost
all the snowdrift formed on the windward side of the two-dimensional and three-dimensional fences, although thesnowdrift also formed along the split streaks on the leeward side of the three-dimensional fence. Our results suggested that the fence should be as long as possible to avoid snowdrifts on roads.
Keywords: Drifting snow, Blowing snow, Snow fence, Lattice Boltzmann method
1 Introduction
Snowdrifts are patchy accumulations of snow resulting
mainly from the redistribution of snow particles on the
ground by drifting snow, which is the horizontal movement of snow particles by creep and saltation on the surface. Snowdrifts can severely affect human activities; forexample, snowdrifts in mountainous terrain can trigger
avalanches and affect the mass balance of water (Lehning
et al. 2006, 2008; Mott et al. 2010; Vionnet et al. 2017).
Moreover, snowdrifts cause traffic disruption because
vehicles tend to stack and stagnate on roads with a snow
depth of 15 cm (Kaneko et al. 2011). Snow fences are oneof the solutions to mitigate the problems caused by snowdrifts, especially on roads in snowy, windy locations.
Snowdrift distribution around snow fences depends
strongly on their design parameters, such as height,*Correspondence: seika@sci.hokudai.ac.jp
1
Graduate School of Science, Hokkaido University, N10W8,
Sapporo 0600810, Japan
Full list of author information is available at the end of the article

thickness, bottom gap, space between panels, penetration
rate of porous fences, and distance from the road (Tabler
1986; Uematsu et al. 1991; Alhajraf 2004). Some studies have examined the onset of drifting snow and snow
accumulation in wind tunnel experiments (Delpech et al.
1998; Okaze et al. 2012; Zhou et al. 2014). Other studies have reported the effect of snow fence design on thesize of snowdrifts on roads based on field observations(Tabler 1980, 1994; Takeuchi 1980; Takeuchi et al. 1986).
However, field observations provide limited opportunities to sample drifting snow events at a certain site, andwind tunnel experiments do not always correspond toreal situations owing to scaling effects. Moreover, it isgenerally costly to obtain sufficient data to design snow
fences using these methods. However, the advancementof computer technology has enabled numerical simulations of snowdrift development to be performed tosearch for an optimal snow-fence design at a target site.
Numerical simulations of drifting snow require computational fluid dynamics (CFD) and the estimation of© The Author(s) 2021. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://creativecommons.org/licenses/by/4.0/.

Tanji et al. Prog Earth Planet Sci(2021) 8:57

snowdrift distribution. Numerical simulations of drifting snow were pioneered in the 1990s by Uematsu
et al. (1991) and Liston et al. (1993). These studies
used a wind simulation based on the Reynolds-averaged Navier–Stokes equations model with turbulence
parameterizations, such as K-theory and the k -ε model.
These models reproduced the snowdrift distribution
around a simple snow fence. Some studies extended
these models to include drifting snow processes due
to multiple snow events that were more complicated(Beyers et al. 2004; Tominaga et al. 2011a, 2011b). The
large-eddy simulation (LES) has been applied to snowtransport simulations to describe turbulence moreaccurately (Okaze et al. 2018; Wang and Jia 2018). For
example, Zwaaftink et al. (2014) combined the LES anda Lagrangian stochastic model and described the temporal and spatial variability of drifting snow with theirmodel. Some studies have considered the momentum
exchange between particles and background wind (Elghobashi 1994). Moreover, turbulent wind is strongly
affected not only by the fixed boundaries, including
topography and artificial obstacles, but also by snow
surfaces that vary temporally due to snowdrifts.
The lattice Boltzmann method (LBM) is a CFD algorithm for quick calculations [see McNamara and Zanetti(1988) for an introduction and Benzi et al. (1992), Qianet al. (1995), and Chen and Doolen (1998) for comprehensive reviews]. In the LBM, the Navier–Stokes equations are replaced with a distribution function equation
called the lattice Boltzmann equation that treats thefluid flow as microscopic fictitious particles in the spacelattice (Chen et al. 1992; Qian et al. 1992). The lattice
Boltzmann equation can be numerically solved by thetranslation of the distribution function and the relaxation
to the equilibrium state. The LBM algorithm is characterized by simpler implementation and higher efficiency inparallel computation than the conventional CFD algorithm (Chen and Doolen 1998). Han et al. (2019) demonstrated that the LBM was three-fold faster than the finite
volume method in 16-core parallel processing. The LBM
has been applied to various fields, such as wind flow inthe urban environment with 1 m resolution (Onodera
et al. 2013), canopy turbulence in neutrally stratified
conditions (Watanabe et al. 2020), flow in porous media(Liu et al. 2016a; b), and flow in blood vessels (Zhang
et al. 2008; Bernasch et al. 2009). In cryology, Wang et al.
(2006) simulated dynamic snowing scenes for various
weather conditions and snow crystal types with LBM.
Lu et al. (2009) also used the LBM to reproduce dendricsnow crystal growth in clouds. These studies suggest that
the LBM is appropriate for modeling blowing snow andsnowdrift distribution; however, no such studies havebeen performed.

Page 2 of 16

In this paper, we develop an LBM model for snowdrift distribution around an artificial snow fence. The
model consists of the CFD module based on the LBM
and a module for calculating the snow particles’ motion
and accumulation following Nishimura and Hunt (2000).
The present work focuses on checking the feasibility ofapplying the LBM to blowing snow and snowdrift modeling in typical experiments. Because the computational
efficiency of the LBM has been demonstrated (Kinget al. 2017; Han et al. 2019), we did not compare the
LBM with other CFD modeling methods. The remaining part of this paper is organized as follows: in Sects. 2
and 3, we describe the model and experiments in detail;in Sect. 4, we show the results of the model simulation,compare them with previous observation and numericalsimulation studies, and discuss the effect of snow fenceson snowdrift distribution; and in Sect. 5 we conclude thepaper.

2 Model2.1 CFD module

The CFD module in the model we developed was basedon a three-dimensional LBM model that delivered microscopic fictitious particles to 19 neighbor nodes, usuallyreferred to as a D3Q19 configuration with the singlerelaxation time collision operator (Fig. 1). This LBM configuration has been established in a sufficiently accurate
turbulent simulation (Onodera et al. 2013; King et al.
2017; Deiterding and Wood 2016; Noh 2019), althoughother sophisticated methods have been proposed toreduce the error (Geier et al. 2009, 2015; Suga et al. 2015).
The discretized lattice Boltzmann equation is a prognostic equation of the distribution function for particles in
the i th direction, fi,

Fig. 1 Schematic of the particle velocity vector, c i , in the D3Q19
configuration of the lattice Boltzmann modelTanji et al. Prog Earth Planet Sci(2021) 8:57

fi (r + c i �t, t + �t) = fi (r, t) + �i �t,

Page 3 of 16(1)

where r is the particle position vector, t is the timeincrement, i is the collision term, and c i is the particle
velocity vector.
The Bhatnagar–Gross–Krook approximation reduced
collision term i to the relaxation to the equilibrium stateof the distribution function (Chen et al. 1992; Qian et al.
1992),
1
�i = − fi − 
fi .
(2)
τ
Here, τ is the relaxation time as a function of viscosity ν∗,

τ=

3ν∗
1
+
2 c2 �t(3)

where c is the discrete speed, defined as x
t , with spatial
increment x. ν∗ is defined as ν∗ = ν0n + νt , where ν0n is
the non-dimensional viscosity of the air (4.0 × 10−6) and νtis the eddy viscosity, given in Eq. (7). 
fi is the equilibrium
distribution function for particles in the i th direction, calculated in the D3Q19 configuration as|u|2(c i • u) (c i • u)2
fi = wi ρ⌈1 +
+
−
⌉,
cs2
2cs4
2cs2

(4)

where




1
3 (i = 0)1
wi =(i = 1, 2, . . . , 6) , 118
(i36 = 7, 8, . . . , 18)(5)√
the lattice speed of sound cs = 1/ 3 , and ρ and u are themacroscopic density and velocity, respectively, calculated
in the D3Q19 configuration as

ρ=

18

i=0

18

fi , andu =1
fi c i .
ρ

(6)

i=0

The dimensional variables are transferred from thenon-dimensional distribution function by multiplying thenon-dimensional value by 50.0 to give the same Reynold
numbers.
The sub-grid scale parameterization (Feng et al. 2007;
Onodera et al. 2013; Wang et al. 2014; Suga et al. 2015) was
implemented to estimate the eddy kinematic viscosity, νt .
We used the Smagorinsky model (Smagorinsky 1963), in
which νt is related to velocity gradient tensor S by 
νt = C�2 S,
(7)

where Smagorinsky coefficient C = 0.12 (Tominaga et al.
2008; Okaze et al. 2021) and C = 60 in the damping zone,
which consists of 15 grids from the outlet boundary todamp the numerical oscillation near the outlet (Inagaki
etal. 2017).  is the cubic root of the local mesh volume.
S was estimated in the D3Q19 configuration as 3 3 18 

 
 
3
2
S  =
fi cij cik ,fi − 
(8)2ρτ
k=1 j=1 i=0

where cij and cik are the j th and k th components of c i,
respectively.
The model domain was a finite channel in the threedimensional space spanned by the wind direction, x , thehorizontal direction normal to x (or fence direction), y ,and the vertical direction, z . Hereafter, for convenience,
the negative and positive ends of x in the domain are
called the western and eastern boundaries, those of y arecalled the southern and northern boundaries, and thoseof z are called the bottom and top boundaries, respectively. The horizontal direction is rotationally invariant inthis system.
The LBM represents the western boundary  condition with inflow u0 (x, y, z) = (u0 y, z , v0 y, z , w0 (y, z))(Sect. 2.2) asρu0,
3
ρ(u0 + v0 ) ρv0
f4 + f17 + f18 − f3 − f15 − f16++,
f7 = f10 +2
63
ρ(u0 − v0 ) ρv0
f3 + f15 + f16 − f4 − f17 − f18
f8 = f9 +
+
−
,
2
6
3
ρ(u0 + w0 ) ρw0
f6 + f16 + f18 − f5 − f15 − f17++
f11 = f14 +
, and
2
6
3
ρ(u0 − w0 ) ρw0
f5 + f15 + f17 − f6 − f16 − f18+−.
f13 = f12 +263
f1 = f2 +(9)Tanji et al. Prog Earth Planet Sci(2021) 8:57

Page 4 of 16

On the other side, the LBM represents the free-flow
condition at the eastern boundary (Hecht and Harting2010) asue,
3ue
f4 + f17 + f18 − f3 − f15 − f16− ,f10 = f7 −46ue
f3 + f15 + f16 − f4 − f17 − f18− ,f9 = f8 −46ue
f6 + f16 + f18 − f5 − f15 − f17− , andf14 = f11 −46
f5 + f15 + f17 − f6 − f16 − f18ue− ,f12 = f13 −
2
6
(10)
where
f2 = f1 −ue = − 1 + f0 + f3 + f4 + f5 + f6 + f15 + f16
+ f17 + f18 + 2 f1 + f7 + f8 + f11 + f13 .

(11)

The southern and northern boundaries are periodic. Distribution functions f3 , f7, f9, f15, and f16 at thesouthern boundary are equal to those at the northernboundary; and distribution functions f4 , f8, f10, f17, andf18 at the northern boundary are equal to those at thesouthern boundary. The top boundary was a free-slip
boundary (i.e., dudz = 0). The LBM representation was

f6 = f5 , f13 = f11 , f14 = f12 , f16 = f15 , and f18 = f17 .
(12)
The bottom boundary was no-slip (i.e., u = 0), so that
the LBM represented it as bounce-back,f5 = f6 , f11 = f14 , f12 = f13 , f15 = f18 , and f17 = f16 .
(13)
The boundary on the fence was also no-slip as bounceback and was written similarly to Eq. (13).
2.2 Inflow turbulence generation

We generated the artificial inflow turbulence andimposed it as the inflow on the western boundary. The
inflow turbulence was generated as two-dimensional
digital-filtered random data by controlling the time andspatial autocorrelations of the resultant inflow turbulence (Okaze and Mochida 2017; Xie and Castro 2008).
The total inflow, u0 (y, z, t), from the western boundary
was divided into the time average, u0 , and the deviation′from the time average, u0. We assumed that the model
domain was the constant flux layer. The wind directionof the time-averaged vector was only eastward, and thewind speed followed the logarithmic profile of
 
z
u∗
,
ln
�u0 (z)� =
(14)
κ
z0
where u∗ is the friction velocity (a parameter to be given),z0 is the roughness length for flat snow surface (0.1 mm;
Hiroo and Ishida 1973), and κ is von Karman’s constant(0.4). Using other common assumptions (Okaze andMochida 2018), Reynolds stress tensor R for the inflow
turbulence was parameterized as


′ ′
′ ′
′
′ 
�u0 u0 � �u0 v0 � �u0 w0 �
10/3 0 −1
′ ′
′ ′
′
′
2
R =  �v0 u0 � �v0 v0 � �v0 w0 �  = u∗  0 5/3 0 .
′ ′
′ ′
′
′
−1 0 5/3
�w0 u0 � �w0 v0 � �w0 w0 �(15)
The deviation
average was computed as ∼ from the time ∼′ 
u0 y, z, t =R (y, z, t), where R is the Cholesky decomposition of R and  is the numerical solution of the stochastic equation (Xie and Castro 2008) ofFig. 2 Schematic of the calculation domain and initial wind vectors
in the experiment with a no fence, b a two-dimensional fence, and ca three-dimensional fence







�t
� y, z, t + �t = � y, z, t exp −
T


1/2


2�t
+ ψ y, z, t + �t 1 − exp −
,
T

(16)

Tanji et al. Prog Earth Planet Sci(2021) 8:57

Page 5 of 16from the initial condition of � y, z, 0 = ψ y, z, 0 . Here,
T is the characteristic timescale and ψ is the digitalfiltered normal random numbers that satisfy a spatial
autocorrelation with characteristic length L (Klein et al.
2003).ψ y, z, t was generated at each time step with a different random number. Following the modified Prandtl
theory (Okuma et al. 1996), L was parameterized as

friction velocity was high enough to lift them (Shaoand Li 1999; Nemoto and Nishimura 2004). Based on
this result, we assumed that a snow particle fell on thebottommost level and occupied the first grid-cell when
the friction velocity on the snow surface was below thethreshold. The friction velocity, u∗ (m s­ −1), on the snow
surface was estimated with a wall function by a twolayer model in Werner and Wenglem (1991) as�2
2ν0 |u(zb )|ν0
for |u(zb )| ≤ 2z
A 1−B

zb
b
� 1
u∗ = �
� �B� �1+B
1+B
1+B
2

ν0
1+B ν0 1−B A 1−B ν0|u(zb )|
, for |u(zb )| ≥ 2z
A 1−B
+ A z
2
z
b

L=

1
κz.
3

b

b

(17)

Then, by assuming Taylor’s hypothesis of frozen turbulence, characteristic timescale T was parameterized
as

T=

1
κz�u0 (z)�.
3(18)

The inflow was numerically generated with the same grid
spacing and the same time interval as the CFD module.
2.3 Snow particle module

The snow particle module in the model followed
Nishimura and Hunt (2000) and Nemoto andNishimura (2004). Assuming that drifting snow particles were spherical, made of ice, electrically neutral,and not driven by the lift force, the equation of motion
for the particles is written as




dup
3 ρa
=−
Cd VR up − u − gk,(19)dt
4 ρp dwhere up (m ­s−1) is the particle
velocity
vector, u (m ­s−1)is the wind vector, VR = up − u , g is gravity (9.8 m ­s−2),
ρa and ρp are the densities of air (1.34 kg ­m−3) and theparticle (910 kg ­m−3), respectively, and d is the particle
diameter (100 µm; Nishimura et al. 2014). Cd is the dragcoefficient for the particle (White 1974), calculated as

Cd =

6
24ν0
+
+ 0.4,VR d
1 + VR d/ν0

(21)

(20)

where ν0 is the viscosity of the air (10−5 m2 s−1). The
terminal fall velocity of a snow particle was estimatedfrom the vertical component of Eqs. (19) and (20) as
0.30 m ­s−1.
Observation results indicated that accumulated snow
particles jumped out of the snow surface when thewhere zb is the height of the bottommost grid just above
the snow surface, A = 8.3, and B = 1/7. This is a different definition of the friction velocity from the inflow generation. The wind velocity at the first grid point variedin time, and then the wall unit in the first grid point wasalso changed. This approach automatically considered
the linear and power law distributions with an instantaneous wall unit calculated with the wind speed in thefirst grid point. The threshold of the friction velocity, u∗t
was computed following Bagnold (1941) and Clifton et al.
(2006) asρp − ρau∗t = 0.2
gd ∼ 0.163 m s−1 .
(22)ρa
We assumed that no snow particles accumulated when
the friction velocity was larger than this threshold. The
aerodynamic entrainment, rebound, splash, mass loss by
sublimation and disruption, the drag force on the fluid,and coalescence of snow particles by the collision were allneglected.

3 Experiments
We performed three experiments: with no fence, with a
two-dimensional fence, and with a three-dimensional
fence of 1.5 m in width (Fig. 2). The fence was set 4 mfrom the western boundary and was centered in thechannel. The fence was non-porous, solid, and the thickness was 0.1 m and the height was 1 m. The following
model settings were independent of the presence of thefence. The channel size was 15.75 × 5 × 5 m, with the grid
spacing at 0.05 m. The origin of the x-coordinate was 4 mfrom the western boundary, at the position of the fence.
In the CFD calculation, the integration time was 30 s
and the time interval was 1 ms; the results were sampled
every 0.02 s from the initial time. In the generation of the
inflow turbulence, the 40 s data was stored with an interval of 4 ms. The inflow turbulence was generated withthe friction velocity, such that the time-averaged windTanji et al. Prog Earth Planet Sci(2021) 8:57

Page 6 of 16

Fig. 3 Vertical profiles of a snow concentration and b snow volume flux of the initial condition on the western boundary plane of the snow particle
modulespeed was 6 m ­s−1 westerly at 10 m aloft in the log profileand the initial wind profile over the calculation domain
was imposed as the same value. We used 30 s data forthe generated inflow turbulence after the spin-up as thewestern boundary condition after linear interpolation as
every 1 ms.
We numerically integrated the equation of motionof snow particles (Eq. 19) with a time interval of 1 ms.
The particles were uniformly distributed on the westernboundary plane per 5.0 cm horizontally and per 2.5 cm
vertically at the initial time, but we assumed that a single
particle in the snow module represented snow mass corresponding to snow volume flux depending on the height(Fig. 3b) as

vf (z) = αn(z)�u0 (z)�
.
ρp

(23)

Here, α = 1500 to accelerate the snow accumulation,and snow concentration n(z) (g m­ −3) is given (Fig. 3a;
Shiotani 1953; Matsuzawa and Takeuchi 2002) by z − 0.30 κu∗,n(z) = min 30, 30
(24)
0.15
where the friction velocity u∗ was estimated by Eq. (14).
The snow particles’ motion was driven by the sampled
time segments of wind data obtained from each experiment with or without a fence; several snow particle integrations were conducted with the CFD output startingfrom a different time. The initial wind profile based on the
CFD module was given to the snow module every 0.1 s
and the wind profile in the snow module was renewed
every 0.02 s. The integration ended when all snow particles had fallen or flowed out of the calculation domain.
Accumulated snow particles gave a height on the gridFig. 4 a Time series of wind speed for the no-fence experiment at
the center of the calculation domain from 0 to 30 s. b Vertical windspeed profile for the no fence experiment at (x, y, t) = (3.5 m, 2.5 m,
0 s) (dotted line) and at (x, y, t) = (3.5 m, 2.5 m, 30 s) (solid line)

corresponding to the snow volume per grid area. Snow
particles within the snowdrift did not move, collapse,
melt, or sublimate. The estimated snowdrift was spatiallysmoothed by averaging data with their neighbors.

Tanji et al. Prog Earth Planet Sci(2021) 8:57

Page 7 of 16

4.2 No‑fence experimentFig. 5 Vertical profile of six components of Reynolds stress a in theartificial inflow turbulence at the center of the y-axis and b in the′ ′no-fence experiment at (x, y) = (3.5 m, 2.5 m). Solid lines are �u u �
′ ′
′ ′
′ ′
(black), �v v � (blue), and �w w � (pink). Broken lines are �u w � (black),′ ′′ ′�u v � (blue), and �v w � (pink). The grey dotted line indicates the
′ ′
′ ′
′ ′
′ ′
′ ′
′ ′
target values of �u w �, �u v � and �v w �, �v v � and �w w �, and �u u �from left to right

4 Results
4.1 Artificial inflow turbulenceFirst, we checked whether the inflow turbulence followed
the target wind profile (Eq. 14) and its corresponding
Reynolds stress tensor (Eq. 15). The time-averaged inflow
turbulence had a vertical profile quite similar to the target log profile with 6 m ­s−1 westerly at 10 m height (notshown). The Reynolds stresses also agreed with the pre′ ′′ ′′ ′scribed values of �u u � ∼ 10/3u2∗,�v v � ∼ �w w � ∼ 5/3u2∗,
′ ′
′
′
′
′
�u w � ∼ −u2∗, and �v u � ∼ �v w � ∼ 0 (Fig. 5a).

A reference experiment without a fence was performed
to demonstrate the validity of the model and to check the
spin-up time. Figure 4a displays the time series of thewind speed at the center of the calculation domain overthe integration period. The calculation was computationally stable during the period. Moreover, because the initial profile occupied the whole calculation domain, it
almost took 2.6 s to run the information from the western boundary to the eastern end. This time span was the
spin-up period of the model. The effect of the inflow turbulence was observed at the central position after half ofthe period (Fig. 4a). The wind data calculated in the CFD
module was discarded before 10 s, and the data segment
selected after 10 s was input into the snow module.
Therefore, we used 201 ensemble members to calculate asnowdrift amount and the snow particles’ trajectories.
The vertical profile of the wind speed at the final time
step is shown in Fig. 4b. The instantaneous westerly windfluctuated around the logarithmic low imposed as themean value at the inflow boundary. Figure 5 shows thesecond-order turbulent statistics at the inflow boundary
in Fig. 5a and the center of the calculation domain′ ′(x = 3.5 m, y = 2.5 m) in Fig. 5b. Although �u u � decreased′ ′slightly and �w w � increased above z = 1.0 m at the centerof the domain, the turbulent
 kinetic energy, defined as
′ ′
′ ′
′ ′
1
2 �u u � + �v v � + �w w � , was comparable with that at
the inflow boundary. The distribution of the three components of the normal stress assumed at the inflow
boundary was only slightly modified. Due to the effect of′ ′the no-slip condition at the ground surface, �w w � wasFig. 6 Snowdrift profile for the no-fence experiment in the cross section along y = 2.5 m, as the sum of the amount of snow accumulation for eachtrajectory calculation

Tanji et al. Prog Earth Planet Sci(2021) 8:57

Page 8 of 16

Fig. 7 Snapshots of the wind vectors (vectors) and vorticity (shading) around the two-dimensional fence in the cross section along y = 2.5 m at a
10 and b 26 s. The solid line at x = 0 shows the two-dimensional fence′′

decreased near the surface, which could decrease �u w �,
especially near the surface at the center of the domain.
′ ′However, �u w � at the center was still half the prescribedvalue at the inflow boundary.
Figure 6 shows the snowdrift distribution as the sumof the amount of snow accumulation in each trajectory
calculation. According to Eqs. (19) and (24), the terminal
fall velocity of snow particles was 0.3 m ­s−1, and snow
particles below a height of 0.5 m at the initial position
had a large volume flux. Snow particles were advectedby the background flow of ~ 4 m ­s−1 until deposition.
Hence, most of the snow particles fell within 7 m of thewestern boundary (Fig. 6). The horizontal distributionof the snowdrift depended on the vertical distribution of
the snow volume flux. The height of the snowdrift was
expected to be about 0.2 m in the calculation domain if
the snow volume flux exceeded 6.0 × 10−7 m3 m−2 s−1,
corresponding to the flux at a height of 0.6 m. No moresnow accumulation was possible after the friction velocity exceeded the threshold (Eq. 22). The results indicatedthat we could test the blocking effect of a fence at x = 0
because a snowdrift was formed beyond the point without the fence.

4.3 Two‑dimensional fence experiment

The horizontal wind observed in the no-fence experiment was distorted by the two-dimensional fence normal to the dominant wind direction, which created an
ascending motion on the windward side and a swirling
eddy in the cross section on the leeward side (two snapshots in Fig. 7). The upper-level wind was intensified
around the top of the fence. These wind profile features
were consistent with previous studies (Uematsu et al.
1991; Liu et al. 2016a, b). However, the eddies were generated successively from the leeward side of the fence andflowed downstream following the dominant wind flow,
and the size and position of the eddies in the cross section varied irregularly over time (Fig. 7). At 10 s, a single
eddy extended 1 m from the fence and there was anotherlarger eddy from x = 6 to 8 m (Fig. 7a). The large eddy
flowed in the eastern direction and left the domain after
11 s (not shown). At 26 s, there were large eddies just to
the east of the fence and further downstream from x = 3
to 6 m (Fig. 7b). These larger eddies had a strong reverse
flow near the surface because the intensified wind at the
top of the fence was separated behind the fence. Eddies

Tanji et al. Prog Earth Planet Sci(2021) 8:57

Page 9 of 16

Fig. 8 a Snowdrift profile around the two-dimensional fence in the cross section along y = 2.5 m, as the sum of the amount of snow accumulationof each trajectory calculation. b–d Ensemble trajectories of snow particles at the center of the y-axis around the two-dimensional fence driven by all
the wind data segments. The initial heights of the particles are b 0.3, c 0.4, and d 1.4 m. The position of the fence is shown by a solid line

successively separated from the east side of the fence. The
eddies caused by the fence augmented the fluctuations
with a larger scale and a lager timescale.

The snowdrift distribution in the two-dimensional
fence experiment (Fig. 8a) was different from that in the
no-fence experiment (Fig. 6). Almost all of the snowdrift

Tanji et al. Prog Earth Planet Sci(2021) 8:57

was distributed on the windward side of the fence
because the surface wind speed was attenuated just to thewest of the fence. The highest snowdrift was about 0.5 m
high at x = −1.3m. These results were consistent withprevious studies using a conventional numerical simulation (Uematsu et al. 1991; Alhajraf 2004) and observations (Tabler 1994). On the leeward side of the fence, few
snowdrifts formed except for from x = 5 to 7 m. This distribution was consistent with the first snowdrift development regime in Tabler (1994), which is explained in detailin Sect. 5. There were three reasons for this result. First,
most of the snow particles that went over the fence did
not accumulate just east of the fence because the strong
wind above the fence accelerated the snow particles. Figure 8b–d show the trajectories of three initial positions of
the snow particles driven by all the wind profile data segments. Most of the particles starting from higher heights
went over the fence and were blown through the calculation domain, such as the particles that started from aheight of 1.4 m (Fig. 8d). However, most of the particles
starting from a height of 0.3 m did not get over the fence(Fig. 8b) because these particles reached the grid on thewindward side of the fence that had a friction velocity below the threshold (Eq. 22). Second, the eddies thatwere successively generated on the leeward side of the

Page 10 of 16

fence disturbed the snowdrift development at a certain
point because these transient eddies were generated byinflow turbulence and the snow accumulation in a single
trajectory calculation depended on the transient eddies.
For example, some particles starting from a height of
0.4 m subsequently followed the streamlines of the swirling eddies and fell to the east of the fence (Fig. 8c), butthese particles fell on different grids in each wind profile
data segment. Third, most of the snow particles accumulating on the leeward side of the fence had a small snow
volume flux, such as the particles starting from a heightof 1.4 m. The snowdrift at the point where these particles fell was less than 5 cm high (Fig. 8a) because theseparticles had a small snow volume flux. In contrast, mostof the snow particles starting at a height of 0.3 m with alarge snow volume flux fell on the surface of the windward side. Thus, the two-dimensional snow fence was aneffective obstacle to snowdrift development on the leeward side compared with the snowdrift distribution in
the no-fence experiment (Fig. 6).
4.4 Three‑dimensional fence experiment

The three-dimensional fence greatly changed the windflow around the fence (Figs. 9, 10). In contrast to the twodimensional fence experiment, the wind flow vector wasFig. 9 Snapshots of wind vectors (vectors) and vorticity (shading) around the three-dimensional fence in the cross section along y = 2.5 m at a 10
and b 26 s. The position of the fence is shown by a solid lineTanji et al. Prog Earth Planet Sci(2021) 8:57

Page 11 of 16

Fig. 10 Snapshots of wind vectors (vectors) and vorticity (shading) around the three-dimensional fence on the horizontal plane along z = 0.5 m at
a 10 and b 26 s. The position of the fence is shown by a solid line

Fig. 11 a Snowdrift height around the three-dimensional fence in the cross section along b y = 2.5 m and c x = 2.5 m. The snowdrift is the sum of
the amount of snow accumulation of each trajectory calculation. The position of the fence is shown by a solid line in a and b

Tanji et al. Prog Earth Planet Sci(2021) 8:57

Page 12 of 16

Fig. 13 Snowdrift potential a for no fence, b around thetwo-dimensional fence, and c around the three-dimensional fence.
The position of the fence is shown by a solid line
Fig. 12 Ensemble trajectories of snow particles at a, b the centerof the y-axis and c, d along z = 0.5 m around the three-dimensional
fence driven by all the wind data segments. The initial positions of theparticles are (y, z) = a (2.5 m, 0.4 m), b (2.5 m, 1.4 m), c (1.5 m, 0.5 m),
and d (2.5 m, 0.5 m). The position of the fence is shown by a solid line

three-dimensional in the three-dimensional fence experiment. The reverse flow along the center of the fence was
weaker than that in the two-dimensional fence experiment around the surface on the leeward side aroundx = 2 m at 10 s (Fig. 9a). The low-level wind flow washorizontally diffluent and went round the fence, andthere was a dipole pattern on the leeward side, with aweak wake flow toward the center of the fence on thehorizontal plane at z = 0.5 m (Fig. 10a). The wind flowat 26 s had a weak eddy from x = 1 to 3 m on the vertical plane (Fig. 9b). These features were consistent with aprevious study that simulated wind flow around a threedimensional obstacle with the LBM (Han et al. 2021). On
the horizontal plane, there was still a dipole pattern, but
the turbulent flow 4 m east of the fence was more intense
after a while due to the spread of turbulence generated
by eddies on the horizontal plane on the leeward side
(Fig. 10b).
In the three-dimensional fence experiment, the snowdrift had a two-dimensional distribution (Fig. 11a). The
snowdrift along the center of the y-axis on the windwardside was similar to that in the two-dimensional fence

experiment (Fig. 11b). On the leeward side, snowdrift
was not formed behind the fence along the center of they-axis, but it was formed in the no-fence zone (Fig. 11c).
This snowdrift developed ahead of the split flow because
there were weak flow lines around the surface (Fig. 10).
This snowdrift profile associated with the split flow wasconsistent with previous studies of snowdrifts aroundcubes (Beyers et al. 2004; Okaze et al. 2013). The snow
particles clarified the three-dimensional trajectories forsnow deposition. All of the particles starting from heights
of 0.3 (not shown) and 0.4 m (Fig. 12a) fell on the windward side of the fence, whereas some particles from thesame heights went over the fence in the two-dimensional
fence experiment. This difference was probably caused by
the lower wind speed just above the three-dimensional
fence due to the energy loss accompanied by the horizontal diffluent flow (Fig. 10). The snow particles’ trajectories
on a horizontal plane indicated the snowdrift formation
process along the split streaks on the leeward side. For
example, some snow particles from (y, z) = (2.5 m, 0.5 m)
did not flow over the fence because they collided withthe fence, whereas most of the snow particles from(y, z) = (1.5 m, 0.5 m) did not collide with the fence andfell along the stream of the horizontal diffluent flow
(Fig. 12c, d). These low-level particles formed the forking
snowdrift on the leeward side.

Tanji et al. Prog Earth Planet Sci(2021) 8:57

4.5 Snowdrift potential

We introduced the snowdrift potential, based on 201
pieces of the snowdrift simulation using the different initial background flow sampled from the CFD calculation from 10 to 30 s. For example, a 50% snowdrift
potential meant that about 100 initial wind profiles satisfied the condition of the development of a snowdriftwith a height of > 5 cm, corresponding to a single verticalgrid interval, when an infinite number of snow particles
flowed in from the western boundary. Figure 13 shows
the snowdrift potential for the three experiments. Snowdrift potential in the no-fence experiment was more than
80% over the calculation area because the friction velocity under a height of about 0.3 m was below the threshold. This result indicated that the snowdrift could form
anywhere in the calculation domain from accumulatingsnow particles blowing from the western boundary. In
contrast, snowdrift potential around the two-dimensional fence decreased greatly on the leeward side of thefence, to less than 40% just 4 m to the right of the fence.
The snowdrift potential was low just behind the threedimensional fence but high elsewhere. The higher probability region extended along two arcs from the bordersof the fence to the downstream area. In addition, mostof the snowdrift potential over 4 m behind the fence wasmore than 80%. These results showed that snow fences
could prevent snowdrift formation on the leeward sideof the fences even if a large number of snow particles
were blown by a strong drifting snow event. However, thetwo-dimensional fence was more effective than the threedimensional fence. The three-dimensional fence used inthis study was not sufficiently long, with a width of 1.5 m
compared with a height of 1 m. Therefore, the separation of eddies near the sides of the fence induced a flow
behind the fence, which led to a more region of low windspeed behind the fence compared with that for the twodimensional fence. As a consequence, the small snowdrift potential region for the three-dimensional fence wasshorter than that of the two-dimensional fence, but theregion with a potential of less than 20% extended up to
3 m behind the fence.

5 Conclusions and discussion
We developed the CFD and snow particle modules toevaluate the snow accumulation around snow fences. The
snow particles were driven by the wind flow in the channel sampled from the CFD model experiment with the
LBM. The snow particles’ motion was modeled following
Nishimura and Hunt (2000) and Nemoto and Nishimura(2004) and the accumulation was computed as a functionof the friction velocity in the viscosity layer. We designed
an experiment with no fence, and experiments with a
two-dimensional fence or a three-dimensional fence,

Page 13 of 16

both installed normal to the dominant wind direction.
The snowdrift distributions in both the two-dimensional
fence and the three-dimensional fence experiments werehigh on the windward side of the fence because mostsnow particles from lower levels did not go over the fence
due to the weak wind, and a varying, strong swirling eddyon the east of the fence often blew snow particles fromhigher levels to out of the domain. However, around thethree-dimensional fence, the snow accumulation was
found on the leeward side along the horizontally diffluent
wind flow. The CFD integration in the three-dimensional
fence experiment increased the vorticity activity of thevertical component, whereas that in the two-dimensional
fence experiment increased the vorticity activity of thehorizontal component. The three-dimensional fence
experiment was designed with the fence blocking 30%of the channel width, but the ratio of the fence length
to channel width probably affected which horizontal orvertical vorticity was dominant in the downward energy
cascade.
The snowdrift model in this paper has two main advantages compared with previous models. First, the ensemble simulation in our model is useful for evaluating theeffect of the fence; it enables us to estimate not only thequantitative snowdrift distribution (Figs. 8a, 11), but also
the probability of the snowdrift development (Fig. 13).
This kind of information helps us to warn drivers by capturing low-probability snowdrift formation events thatcan cause traffic disruption. Second, temporal variationof the boundary conditions of the CFD module due tosnowdrift formation can be added easily to our model.
Most previous studies did not update the snow surface
boundary in the calculation because of the high computational cost of conventional CFD algorithms and thetechnical complexity of updating boundary conditions
during the simulation.
Even though this paper focused on model development, we discuss the feasibility of our model on a qualitative reproduction by our model of three regimes ofsnowdrift formation around a solid fence described by
Tabler (1994). In Regime I, the snowdrift developed on
the windward side of the fence with cavity between thefence and the snowdrift, and the snowdrift was 0.6 times
the height of the fence. Our simulation result for the twodimensional fence was consistent with the characteristics of Regime I (Fig. 8a). However, the simulation results
did not reproduce the other regimes of filling the cavity(Regime II) and snowdrift development on the leeward
side of the fence (Regime III). Our model should be parallelized to enable a longer integration time to reproduce
these regimes in a reasonable computational time. Moreover, the boundary condition in the CFD module must be
updated successively, so as to include a possible change ofTanji et al. Prog Earth Planet Sci(2021) 8:57

flow due to a change in snowdrift surface. Furthermore,
the resuspension and redeposition of the particles should
be included in the snow particle module to estimate
snowdrifts more accurately.
To achieve a more accurate snowdrift simulation,other snow motion and accumulation processes must be
included. For example, resuspension of snow particles
from the surface is an important process in drifting snow.
In the present model, this process was implicitly includedin the prohibition of snow accumulation on the surface
by strong wind. However, we did not consider the trajectories and redeposition of resuspended snow particles.
The resuspension processes are aerodynamic entrainment, rebound and splash (Shao and Li 1999; Ammi et al.
2009). Moreover, the initial condition of the snow surface
should be prescribed because snow particles on the surface drift when the friction velocity exceeds the thresholdvelocity.
This paper was limited to experiments with no interactions between snow and wind. In this study, snow
particles were affected by the wind, but the wind was
not affected by the snow particles; thus, the coupling
was one-way. By considering the interaction between
snow and wind as a two-way coupling, the wind velocity is slightly reduced by the momentum exchange
between snow particles and wind. This modification
may change the trajectory of snow particles (Figs. 8b–d, 12) so that they fall short on the windward side of
the fence. Moreover, there is an interaction between
the snow surface and wind flow. The snow accumulation changes the bottom boundary condition in the
CFD calculation. Our model can be extended easily toallow temporal variation of the boundary conditions
because the LBM is simpler than other algorithms andmore suitable for complicated boundary conditions.
Because the snow particles generally accumulated
where the wind speed was low on the windward side

Page 14 of 16

of the three-dimensional fence (Fig. 11), we can easily
presume that the snow surface on the windward side is
asymptotic to a streamline that crosses the top of thefence. We can readily implement this interaction process simply by combining the CFD and the snow particle modules, but this will be addressed in future work.
Although this paper is limited to examining the feasibility of applying the LBM to drifting snow and snowdrift modeling, we can still compare the simulation
results with observations. An Observation of driftingsnow around a long, wide, solid fence in Teshikaga, asmall town in eastern Hokkaido, on February 18, 2019were recorded over several hours (Okaze et al. 2019).
The height and width of the fence were 1 and 6 m,
respectively. The snow depth was measured along aline orthogonal to the center of the fence. The winddirection did not change much and was orthogonal to
the fence. The fence setting was quite similar to that in
the two-dimensional numerical experiment shown in
Sect. 4.3. Comparing the observation with our results(Fig. 14) indicated that the snowdrift distribution in
the numerical experiment and the peak location of thesnowdrift at about x = −1 m were consistent with theobservation However, we cannot compare these results
further because of the lack of high-resolution spatiotemporal observation data. Moreover, a fairer comparison is needed to improve the experiment design
and improve the model components to represent realistically the inflow and surface boundary conditions andtheir interactions with drifting snow particles, such as
rebounding and resuspension.
Abbreviations
CFD: Computational fluid dynamics; LBM: Lattice Boltzmann method; LES:
Large-eddy simulation.
Acknowledgements
We thank Dr. Yousuke Sato for providing insightful comments in the earlystage of this study, and Drs. Masaru Matsuzawa, Yusuke Harada, and Satoshi
Omiya for organizing the blowing-snow observations that encouraged us toperform this study.
Authors’ contributions
ST proposed the topic, conceived and designed the study, developed themain source code of the snowdrift model, analyzed the results, and wrote themanuscript. MI helped in the development of the snowdrift model, interpretation of the results, and writing the manuscript. TO made the artificial inflow
turbulent and helped in the development of the snowdrift model. All authors
discussed the results and read and approved the final manuscript.

Fig. 14 Snowdrift profile around the two-dimensional fence in the
cross section along y = 2.5 m (line) (as shown in Fig. 8a) and thesnowdrift profile observed in Okaze et al. (2019) (blue circles)

Funding
This work is supported by the Environment Research and Technology
Development Fund JPMEERF20192005 of the Environmental Restoration andConservation Agency of Japan, the Integrated Research Program for Advancing Climate Models (TOUGOU) Grand number JPMXD0717935457 from
the Ministry of Education, Culture, Sports, Science, and Technology (MEXT),by JSPS KAKENHI Grants (18K03734, 18H03819, and 19H00963), and by the
Research Field of Hokkaido Weather Forecast and Technology Development(endowed by Hokkaido Weather Technology Center Co., Ltd.).

Tanji et al. Prog Earth Planet Sci(2021) 8:57

Availability of data and materials
The datasets used during the current study are available from the corresponding author on reasonable request.

Declarations
Competing interests
The authors declare that they have no competing interests.
Author details
1
Graduate School of Science, Hokkaido University, N10W8, Sapporo 0600810,
Japan. 2 Faculty of Science, Hokkaido University, N10W8, Sapporo 0600810,
Japan. 3 Center for Natural Hazards Research, Hokkaido University, N9W9,
Sapporo 0608589, Japan. 4 School of Environment and Society, Tokyo Instituteof Technology, 4259 Nagatsuta‑cho, Midori‑ku, Yokohama 2268502, Japan.
Received: 19 March 2021 Accepted: 27 September 2021

References
Alhajraf S (2004) Computational fluid dynamic modeling of drifting particles atporous fences. Environ Model Softw 19:163–170. https://​doi.​org/​10.​1016/​
S1364-​8152(03)​00118-X
Ammi M, Oger L, Beladjine D, Valance A (2009) Three-dimensional analysisof the collision process of a bead on a granular packing. Phys Rev E
79:021305. https://​doi.​org/​10.​1103/​PhysR​evE.​79.​021305
Bagnold RA (1941) The physics of blown sand and desert dunes. Methuen,
LondonBenzi R, Succi S, Vergassola M (1992) The Lattice Boltzmann equation: theory
and applications. Phys Rep 222:145–197. https://​doi.​org/​10.​1016/​0370-​
1573(92)​90090-M
Bernasch M, Fatica M, Melchionna S, Succi S, Kaxiras E (2009) A flexible highperformance lattice Boltzmann GPU code for the simulations of fluid
flows in complex geometries. Concurr Comput Pract Exp 22:1–14. https://​doi.​org/​10.​1002/​cpe.​1466
Beyers JHM, Sundsbo PA, Harms TM (2004) Numerical simulation of threedimensional, transient snow drifting around a cube. J Wind Eng Ind
Aerodyn 92:725–747. https://​doi.​org/​10.​1016/j.​jweia.​2004.​03.​011
Chen S, Doolen GD (1998) Lattice Boltzmann method for fluid flows. Annu Rev
Fluid Mech 30:329–364. https://​doi.​org/​10.​1146/​annur​ev.​fluid.​30.1.​329
Chen H, Chen S, Matthaeus WH (1992) Recovery of the Navier Stokes equations using a lattice-gas Boltzmann method. Phys Rev A 45:5339–5342.
https://​doi.​org/​10.​1103/​PhysR​evA.​45.​R5339
Clifton A, Rüedi JD, Lehning M (2006) Snow saltation threshold measurements
in a drifting snow wind tunnel. J Glaciol 39:585–596. https://​doi.​org/​10.​3189/​17275​65067​81828​430
Deiterding R, Wood SL (2016) Predictive wind turbine simulation with an adaptive lattice Boltzmann method for moving boundaries. J Phys Conf Ser
753:082005. https://​doi.​org/​10.​1088/​1742-​6596/​753/8/​082005
Delpech P, Palier P, Gandemer J (1998) Snowdrifting simulation around Antarctic buildings. J Wind Eng Ind Aerodyn 74–76:567–576. https://​doi.​org/​10.​1016/​S0167-​6105(98)​00051-8
Elghobashi S (1994) On predicting particle-laden turbulent flows. Appl Sci Res
52:309–329. https://​doi.​org/​10.​1007/​BF009​36835
Feng YT, Han K, Owen DRJ (2007) Coupled lattice Boltzmann method anddiscrete element modelling of particle transport in turbulent fluid flows:
computational issues. Int J Numer Methods Eng 72:1111–1134. https://​doi.​org/​10.​1002/​nme.​2114
Geier M, Greiner A, Korvink JG (2009) A factorized central moment lattice
Boltzmann method. Eur Phys J Spec Top 171(1):55–61. https://​doi.​org/​10.​
1140/​epjst/​e2009-​01011-1
Geier M, Schönherr M, Pasquali A, Krafczyk M (2015) The cumulant lattice
Boltzmann equation in three dimensions: theory and validation. Comput
Math Appl 70(4):507–547. https://​doi.​org/​10.​1016/j.​camwa.​2015.​05.​001
Han M, Ooka R, Kikumoto H (2019) Lattice Boltzmann method-based largeeddy simulation of indoor isothermal airflow. Int J Heat Mass Transf
130:700–709. https://​doi.​org/​10.​1016/j.​ijhea​tmass​trans​fer.​2018.​10.​137

Page 15 of 16

Han M, Ooka R, Kikumoto H (2021) Effects of wall function model in lattice
Boltzmann method-based large-eddy simulation on built environment
flows. Build Environ 195:107764. https://​doi.​org/​10.​1016/j.​build​env.​2021.​
107764
Hecht M, Harting J (2010) Implementation of on-site velocity boundary conditions for D3Q19 lattice Boltzmann simulations. J Stat Mech. https://​doi.​
org/​10.​1088/​1742-​5468/​2010/​01/​P01018
Hiroo F, Ishida T (1973) Rate of Turbulent energy dissipation during snow drifting. Low Temp Sci Ser A Phys Sci 31:69–85
Inagaki A, Kanda M, Ahmad HN, 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. Bound Layer Meteor
164:1651–2181. https://​doi.​org/​10.​1007/​s10546-​017-​0249-y
Kaneko M, Matsuzawa M, Watanabe T (2011) Experiments on snowdrift development and vehicle starting limit. Technical research presentation 2011
(in Japanese)
King MF, Khan A, Delbosc N, Gough HL, Halios C, Barlow JF, Noakes CJ (2017)
Modelling urban airflow and natural ventilation using a GPU-basedlattice-Boltzmann method. Build Environ 125:273–284. https://​doi.​org/​10.​1016/j.​build​env.​2017.​08.​048
Klein M, Sadiki A, Janicka J (2003) A digital filter based generation of inflow
data for spatially developing direct numerical or large eddy simulations.
J Comput Phys 186:652–665. https://​doi.​org/​10.​1016/​S0021-​9991(03)​00090-1
Lehning M, Volksch I, Gustafsson D, Nguyen TA, Stahli M, Zappa M (2006)
ALPINE3D: a detailed model of mountain surface processes and its
application to snow hydrology. Hydrol Proc 20:2111–2128. https://​doi.​
org/​10.​1002/​hyp.​6204
Lehning M, Löwe H, Ryser M, Raderschall N (2008) Inhomogeneous precipitation distribution and snow transport in steep terrain. Water Resour Res
44:W07404. https://​doi.​org/​10.​1029/​2007W​R0065​45
Liston GE, Brown RL, Dent JD (1993) A two-dimensional computational modelof turbulent atmospheric surface flows with drifting snow. Ann Glaciol
18:281–286. https://​doi.​org/​10.​3189/​S0260​30550​00116​54
Liu D, Li Y, Wang B, Hu P, Zhang J (2016a) Mechanism and effects of snow
accumulations and controls by lightweight snow fences. J Mod Transp
24:261–269. https://​doi.​org/​10.​1007/​s40534-​016-​0115-5
Liu H, Kang Q, Leonardi CR, Schmieschek S, Narváez A, Jones BD, Williams JR,
Valocchi AJ, Harting J (2016b) Multiphase lattice Boltzmann simulations
for porous media applications. Comput Geosci 20:777–805. https://​doi.​
org/​10.​1007/​s10596-​015-​9542-3
Lu GP, DePaolo DJ, Kang QJ, Zhang DX (2009) Lattice Boltzmann simulation ofsnow crystal growth in clouds. J Geophys Res 114:D07305. https://​doi.​
org/​10.​1029/​2008J​D0110​87
Matsuzawa M, Takeuchi M (2002) A study of methods to estimate visibility
based on weather conditions. Seppyo 64:77–85 (in Japanese)
McNamara GR, Zanetti G (1988) Use of the Boltzmann equation to simulatelattice-gas automata. Phys Rev Lett 61:2332–2335. https://​doi.​org/​10.​1103/​PhysR​evLett.​61.​2332
Mott R, Schirmer M, Bavay M, Grünewald T, Lehning M (2010) Understandingsnow-transport processes shaping the mountain snow-cover. Cryosphere
4:545–559. https://​doi.​org/​10.​5194/​tc-4-​545-​2010
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/​2004J​D0046​57
Nishimura K, Hunt JCR (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/​S0022​11200​00010​14
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 119:9901–
9913. https://​doi.​org/​10.​1002/​2014J​D0216​86
Noh HM (2019) Numerical analysis of aerodynamic noise from pantographin high-speed trains using lattice Boltzmann method. Adv Mech Eng
11(7):1–12. https://​doi.​org/​10.​1177/​16878​14019​863995
Okaze T, Mochida A (2017) Cholesky decomposition–based generation of
artificial inflow turbulence including scalar fluctuation. Comput Fluids
159:23–32. https://​doi.​org/​10.​1016/j.​compf​l uid.​2017.​09.​005
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 Aerodyn 104–106:532–539. https://​doi.​org/​10.​1016/j.​jweia.​2012.​04.​002

Tanji et al. Prog Earth Planet Sci(2021) 8:57

Okaze T, Tominaga Y, Mochida A (2013) Development of new snowdrift model
based on two transport equations of drifting snow density, numericalprediction of snowdrift around buildings using CFD (part 2). J Environ
Eng AIJ 78:149–156 (in Japanese)
Okaze T, Niiya H, Nishimura K (2018) Development of a large-eddy simulation
coupled with Lagrangian snow transport model. J Wind Eng Ind Aerodyn
183:35–43. https://​doi.​org/​10.​1016/j.​jweia.​2018.​09.​027
Okaze T, Kikumoto H, Ono H, Imano M, Ikegaya N, Hasama T, Nakao K, Kishida
T, Tabata Y, Nakajima K, Yoshie R, Tominaga Y (2021) Large-eddy simulation of flow around an isolated building: a step-by-step analysis of
influencing factors on turbulent statistics. Build Environ. https://​doi.​org/​
10.​1016/j.​build​env.​2021.​108021
Okaze T, Mochida A (2018) Generation of artificial Inflow turbulence with digital filter method for LES within built-up environment. In: 25th Symposiumof wind engineering, pp 211–216 (in Japanese)
Okaze T, Niiya H, Omiya S, Sunako S, Nishimura K (2019) Observation of snowdrift around two-dimensional fence during a drifting snow event. In: JSSI
and JSSE joint conference 2019 in Yamagata, p 110 (in Japanese)
Okuma T, Kanda J, Tamura Y (1996) Windproof design of buildings. Kajima
Syuppan-kai (in Japanese)
Onodera N, Aoki T, Shimokawabe T, Kobayashi H (2013) A large-scale LES windsimulation using lattice Boltzmann method for a 10 km × 10 km areain Tokyo. In: High performance computing symposium, pp 123–131 (inJapanese)
Qian Y, d’Humieres D, Lallemand P (1992) Lattice BGK models for Navier-Stokes
equation. Europhys Lett 17:479–484. https://​doi.​org/​10.​1209/​0295-​5075/​
17/6/​001
Qian YH, Succi S, Orszag SA (1995) Recent advances in lattice Boltzmann
computing. Annu Rev Comput Phys 3:195–242. https://​doi.​org/​10.​1142/​
97898​12830​647_​0006
Shao Y, Li A (1999) Numerical modelling of saltation in the atmospheric surface layer. Bound Layer Meteorol 91:199–225. https://​doi.​org/​10.​1023/A:​
10018​16013​475
Shiotani M (1953) Note on the vertical distribution of density of blowing snow
in the snow storm. Seppyo 15:6–9 (in Japanese)
Smagorinsky J (1963) General circulation experiments with the primitive equations. I. The basic experiment. Mon Weather Rev 91:99–164. https://​doi.​
org/​10.​1175/​1520-​0493(1963)​091%​3c0099:​GCEWTP%​3e2.3.​CO;2
Suga K, Kuwata Y, Takashima K, Chikasue R (2015) A D3Q27 multiple-relaxation-time lattice Boltzmann method for turbulent flows. Comput MathAppl 69:518–529. https://​doi.​org/​10.​1016/j.​camwa.​2015.​01.​010
Tabler RD (1980) Geometry and density of drifts formed by snow fences. J
Glaciol 26:405–419. https://​doi.​org/​10.​3189/​S0022​14300​00109​35
Tabler RD (1986) Snow fence handbook (Release 1.0). Tabler and Associates,
Laramie, Wyoming
Tabler RD (1994) Design guidelines for the control of blowing and driftingsnow. Strategic highway research program report SHRP-H-381
Takeuchi M (1980) Vertical profile and horizontal increase of drift-snow transport. J Glaciol 26:481–492. https://​doi.​org/​10.​3189/​S0022​14300​00109​96
Takeuchi M, Ishimoto K, Nohara T, Fukuzawa Y (1986) Dynamic threshold windspeed for suspension. In: Proceedings of the annual meeting of theJapanese Society of Snow and Ice (in Japanese)
Tominaga Y, Mochida A, Murakami S, Sawaki S (2008) Comparison of various
revised k–ε models and LES applied to flow around a high-rise building

Page 16 of 16

model with 1:1:2 shape placed within the surface boundary layer. J WindEng Ind Aerodyn 96:389–411. https://​doi.​org/​10.​1016/j.​jweia.​2008.​01.​004
Tominaga Y, Mochida A, Okaze T, Sato T, Nemoto M, Motoyoshi H, Nakai S, Tsutsumi T, Otsuki M, Uamatsu T, Yoshino H (2011a) Development of a systemfor predicting snow distribution in built-up environments: combining
a mesoscale meteorological model and a CFD model. J Wind Eng Ind
Aerodyn 99:460–468. https://​doi.​org/​10.​1016/j.​jweia.​2010.​12.​004
Tominaga Y, Okaze T, Mochida A (2011b) CFD modeling of snowdrift aroundabuilding: an overview of models and evaluation of a new approach.
Build Environ 46:899–910. https://​doi.​org/​10.​1016/j.​build​env.​2010.​10.​020
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, Lac C, Bouvet F, Guyomarc’h G (2017) Highresolution large eddy simulation of snow accumulation in Alpine terrain.
J Geophys Res 122:11005–11021. https://​doi.​org/​10.​1002/​2017J​D0269​47
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-​3841-​2018
Wang CB, Wang ZY, Xia T, Peng QS (2006) Real-time snowing simulation. Vis
Comput 22:315–323. https://​doi.​org/​10.​1007/​s00371-​006-​0012-8
Wang X, Shangguan YQ, Onodera N, Kobayashi H, Aoki T (2014) Direct numerical simulation and large eddy simulation on a turbulent wall-boundedflow using Lattice Boltzmann Method and multiple GPUs. Math Probl
Eng. https://​doi.​org/​10.​1155/​2014/​742432
Watanabe T, Shimoyama K, Kawashima M, Mizoguchi Y, Inagaki A (2020)
Large-eddy simulation of neutrally-stratified turbulent flow within andabove plant canopy using the central-moments-based lattice Boltzmann
method. Bound Layer Meteorol 176:35–60. https://​doi.​org/​10.​1007/​
s10546-​020-​00519-8
Werner H, Wenglem H (1991) Large-eddy simulation of turbulent flow over
and around a cube in a plate channel. In: Proceedings of 8th symposium
on turbulent shear flows, pp 155–158. https://​doi.​org/​10.​1007/​978-3-​642-​
77674-8_​12.
White F (1974) Viscous fluid flow. McGraw-Hill Inc., New YorkXie ZT, Castro IP (2008) Efficient generation of inflow conditions for large
eddy simulation of street-scale flows. Flow Turbul Combust 81:449–470.
https://​doi.​org/​10.​1007/​s10494-​008-​9151-5
Zhang J, Johnson PC, Popel AS (2008) Red blood cell aggregation and dissociation in shear flows simulated by lattice Boltzmann method. J Biomech
41:47–55. https://​doi.​org/​10.​1016/j.​jbiom​ech.​2007.​07.​020
Zhou X, Hu J, Gu M (2014) Wind tunnel test of snow loads on a stepped
flat roof using different granular materials. Nat Hazards 74:1629–1648.
https://​doi.​org/​10.​1007/​s11069-​014-​1296-z
Zwaaftink CDG, Diebold M, Horender S, Overney J, Lieberherr G, Parlange MB,
Lehning M (2014) Modelling small-scale drifting snow with a Lagrangianstochastic model based on large-eddy simulations. Bound Layer Meteorol 153:117–139. https://​doi.​org/​10.​1007/​s10546-​014-​9934-2

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.