Tanji et al.
Progress in Earth and Planetary Science(2023) 10:67
https://doi.org/10.1186/s40645-023-00599-3

Progress in Earth andPlanetary Science

Open Access

RESEARCH ARTICLE

Comparison of the LBM snowdrift model
output with the observation results
Seika Tanji1* , Masaru Inatsu2,3 and Yusuke Harada4

Abstract
In this work, snowdrift experiments which are equivalent to one drifting snow event are performed by the snowdrift
model. The model consisted of the computational fluid dynamics part of the large-eddy simulation with the lattice
Boltzmann method and the drifting snow part of the conventional advection algorithm for representative Lagrangianparticles. The observed vertical wind profile of a 4 h drifting snow event in Teshikaga Town was used as the inflow
boundary conditions in the model to compare the results of the snowdrift estimated by the model and the observed
snowdrift distribution. Parallelization enabled us to simulate the snowdrift distribution in a realistic domain and on thetime scale of a single drifting snow event. We demonstrated that the upgraded model could quantitatively reproduce
the height and position of the observed snowdrift along the center of a three-dimensional fence.
Keywords Blowing snow, Drifting snow, Lattice Boltzmann method, Snowdrift, Snow fence
1 Introduction
Drifting snow is a phenomenon in which surface snow
particles or snowfall particles are blown by a strong windand often form snowdrifts in places where the windspeed is sufficiently low due to buildings and landforms.
The numerical simulation of snowdrifts has overcome thelimited observation opportunities of drifting snow eventsand the rough assumption of similarity in wind tunnel experiments. The most important issue in snowdrift
simulation is resolving wind profiles at and above thesurface because snow particles suspended in the air are
advected by the background wind. Uematsu et al. (1991)
and Liston et al. (1993) pioneered snowdrift simulation
with wind profiles calculated by the Reynolds-averaged*Correspondence:
Seika Tanjiseika@storm.dpri.kyoto-u.ac.jp
1
Disaster Prevention Research Institute, Kyoto University, Gokasho, Uji,
Kyoto 6110011, Japan
2
Faculty of Science, Hokkaido University, N10W8, Sapporo 0600810,
Japan
3
Center for Natural Hazards Research, Hokkaido University, N9W9,
Sapporo 0608589, Japan
4
Civil Engineering Research Institute for Cold Region, 1‑34 1‑13 Hiragishi,
Tohohira‑Ku, Sapporo 0628602, Japan

Navier–Stokes equations model with turbulence parameterizations. A CFD (computational fluid dynamics)
model for snowdrift simulation was developed into the
LES (large-eddy simulation) with the finite differential
method (Zwaaftink et al. 2014; Okaze et al. 2018; Wangand Jia 2018).
A stumbling block for snowdrift simulation is demanding CFD calculations. CFD calculations on the order ofcentimeters are difficult for reproducing the observed
snowdrift for model validation on a realistic spatial
scale and event time scale. Tanji et al. (2021) developed
a snowdrift model that contained the CFD part of the
LES with the LBM (lattice Boltzmann method; McNamara and Zanetti 1988) and the drifting snow part ofthe conventional advection algorithm for representative
Lagrangian particles. They showed that this model could
be used for simulating snowdrifts based on the reasonable results from their idealized experiments. In addition,
the LBM in this model had a simpler implementation andhigher parallel computation efficiency than the conventional CFD algorithms (Chen and Doolen 1998; Han et al.
2019). Hence, parallelizing the model should enable us toperform snowdrift experiments with a full three-dimensional CFD computation for a sufficiently large domain© The Author(s) 2023. 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. Progress in Earth and Planetary Science(2023) 10:67

and long-time integration, which has not been achieved
yet.
Here, this paper aims to enable to conduct a long-run
numerical simulation which is equivalent to one drifting snow event and reproduce an observed snowdrift
distribution quantitatively by parallelizing the model
developed by Tanji et al. (2021). The snowdrift part ofthe model was also upgraded to allow quantitative treatment of drifting snow particles. We added resuspension
and rebound processes that seemed non-negligible, but
excluded splash (Sugiura and Maeno 2000) and sublimation that likely contributed little to snowdrift formationin the calculation domain and time scale in this study(Okaze et al. 2018; Palm et al. 2017; Zwaaftink et al.
2013). The experiments in this study were limited to thetime scale for snow accumulation upstream of obstacles, here a solid fence, so that we ruled out the interaction between fluid dynamics and snow surface change.
Because there was no turbulent wind observation in
the full domain, we imposed the laminar flow fit to thesingle-point anemometer observations as the upstream
boundary after we checked that the turbulent inflow did
not substantially affect the snowdrift results in advance.
The remaining part of this paper is organized as follows. We introduce the observation data in Sect. 2 andthe model upgraded from Tanji et al. (2021) in Sect. 3; we
describe the experiment settings in Sect. 4; we show the

Page 2 of 12

simulation results of the snowdrift distribution and compare it with the observations in Sects. 5; and we conclude
the paper in Sect. 6.

2 Observation data
We used observation data acquired in Teshikaga, Hokkaido, Japan (43.502° N and 144.467° E; Fig. 1; Okaze et al.
2019). This site was a flat field with no obstacle over 500
m on the windward side of the observation point. On thesite, wind speed and direction were observed every 10
min with three-dimensional ultrasonic anemometers atfour levels of the observation tower, at heights of about 1,
1.5, 3, and 7 m above the surface. The number of driftingsnow particles passing the specific area was also measured at the same levels by SPCs (snow particle counters)with every second output. The solid, non-porous fence
was 6 m long, 4 mm thick, and 1 m high and was set in
the field (Fig. 2a). The sensors of ultrasonic anemometersand SPCs and the fence are located on almost the sameline perpendicularly to the dominant wind direction.
On February 18, 2019, strong wind and drifting snow
were recorded in the morning, and a snowdrift wasformed around the fence. The day before the driftingsnow event, there had been snowfall with light wind,
and the field was uniformly covered with snow about 2
cm depth. However, no snowfall was recorded, while thedrifting snow was developing on February 18 according

Fig. 1 Map of northern Japan and surrounding area. The blue circle denotes the observation site in Teshikaga Town

Tanji et al. Progress in Earth and Planetary Science(2023) 10:67

Page 3 of 12

Fig. 2 a Photograph showing the observation tower loaded with anemometers and snow particle counters and the fence at the observationsite. b Photograph showing the measurement of snowdrift height in the cross section along the center of the fence. Vectors show the main winddirections at the site. These photographs were taken on another day, not related to the event that we examined

to the observations from a heated rain gauge at an AMeDAS (Automated Meteorological Data Acquisition System) point near the site (43.51° N and 144.45° E). On
the day following the drifting snow event, the snowdrift
height was measured in the cross section normal to thefence at its center (Fig. 2b).

3 Method
3.1 CFD module

The CFD module was the same as in Tanji et al. (2021),
but the program was fully parallelized. The module
was an LBM LES model with the D3Q19 configuration,
which calculated the friction velocities of 3-dimensional
19 directions. The single relaxation time was imposed
as a collision operator. The subgrid scale parameterization was applied to estimate the eddy kinematic viscosity (Smagorinsky 1963; Feng et al. 2007; Onodera et al.
2013; Wang et al. 2014; Suga et al. 2015). Different from
Tanji et al. (2021), the non-dimensional viscosity of the
air which was changed to 2.0 × 10−6 to simulate the flow
with a higher Reynolds number.
We set the model domain to a three-dimensional channel, which was almost the same with the three-dimensional fence experiment in Tanji et al. (2021). We defined
wind direction, x axis, as from the western to easternboundaries. The model coordinate is set as x axis for thewind direction, y axis for the fence direction normal to
the wind and z axis for the vertical direction. The boundary conditions were the same as Tanji et al. (2021), but theinitial and the inflow conditions on the western boundary
were replaced with observation data (Sect. 4.1). The east,bottom, and top boundaries were free flow, no-slip, andfree-slip conditions, respectively. The north and south
boundaries conditions were imposed on cyclic.

3.2 Snow particle module

The snow particle module consisted of moving, deposition, resuspension, and rebound; the latter two wereadded to the model developed in Tanji et al. (2021).
The resuspension is the process that snow particles on
the surface are peeled off by the strong wind, and therebound is the process that snow particles bound on thesurface mechanically.
The moving process described the blown snow particles in the air by the wind calculated by the CFD modulewith referred to Nishimura and Hunt (2000). The method
assumed that the drifting snow particles are spherical andcalculated the particles’ speed related to the wind flow
and drag coefficient. We set the particle diameter to 135µm and density of particles to 910 kg ­m−3 (Nishimura andHunt 2000). The given diameter was close to the average
diameter of 135 µm observed at heights of 1 and 1.5 m at
the site because snow particles following at lower levels
contribute to snowdrift formation (Tanji et al. 2021).
The resuspension process was activated when the windspeed near the surface was strong enough to peel off
accumulated snow particles on the surface. The thresholdvalue of the resuspension is decided by the friction velocity following Clifton et al. (2006) asu∗t = 0.2

ρp − ρa
gd,
ρa

(1)

where g is the gravity (9.8 m ­s−2), ρa is the densities ofair (1.34 kg ­m−3), d is the particle diameter, and ρp isthe density of snow particles. The threshold value wasabout 0.189 m ­s−1 in the condition of this study. The
friction velocity on each grid was estimated by Wernerand Wenglem (1991) (see Tanji et al. 2021). A snow particle was launched from a grid where the friction velocity exceeded the threshold value. The initial velocity ofTanji et al. Progress in Earth and Planetary Science(2023) 10:67

the resuspended particle was equal to the wind speed at
one grid above the resuspended point. A single particle
represented snow mass corresponding to the volume ofsnow on the grid (Sect. 4.1). A launched particle was suspended in the air, and the volume was released when theparticle accumulated on the surface.
A snow particle that reached the surface proceeded to
the rebound or the deposition process. The process after
a collision was determined by the kinetic energy of thesnow particle. If the particle velocity after the collision
with the surface was enough for the particle to bounce by
one grid height in the CFD calculation, 5 cm, the particle
rebounded. Otherwise, the particle proceeded to the deposition process and stayed unless the friction velocity there
was strong enough to activate the resuspension process.
Elevation angle θr , azimuth angle ϕr , and restitution coefficient er were estimated when particles collided with thesurface (Fig. 3). Although Okaze et al. (2018) suggested thatthese parameters follow the normal distributions, we gave
the fixed values as a function of incident angle θin as elevation angle θr = 20◦ + 0.19θin, azimuth angle ϕr = 0, andrestitution coefficient er = 0.87 − 0.62sinθin, which were
the averaged value in Okaze et al. (2018). We also assumed
that a particle would consecutively rebound 50 times atmost because a few particles continuously rebounded in
the same grid.

4 Experiments
4.1 Inflow condition

The imposed wind profile as the initial condition and thewestern boundary condition in the CFD module was thelaminar logarithmic profile with the roughness length fora flat snow surface z0 = 0.1(mm) (Nishio and Ishida 1973)
 
z
u∗
,
ln
u0 (z) =(2)
κz0
where von Karman’s constant κ = 0.4 . The friction velocity was determined by observed wind data using thePage 4 of 12

least squares method by the following procedure. First,
we chose the period of the drifting snow event from theobserved wind data and the amount of drifting snow
data and averaged the wind data at each height over theperiod. Second, we decided the friction velocity in Eq. (2)with the minimum root mean square error between the
averaged observation data and the logarithmic profile.
Finally, we obtained the logarithmic profile, which wasclose to the observation results. In this paper, we did not
apply artificial inflow turbulence as the inflow condition
because the inflow turbulence did not substantially affectthe results for snowdrifts on the windward side of thefence, as shown in Additional file 1.
The amount of inflow snow particles from the westernboundary was also estimated by the logarithmic profile
decided above. We did not directly use SPCs data as theamount of inflow snow particles. This was because thevertical number of SPCs was insufficient at lower level
where inflow snow particles largely contributed to thesnowdrift formation (Tanji et al. 2021). Instead, following
Tanji et al. (2021), snow particles were initially arranged
in a grid of 5 cm in horizontal direction by 2.5 cm invertical direction on the inflow plane. A single particle in this model represented snow mass corresponding
to snow volume flux vf (m3 m−2 s−1 ) depending on theheight with referred to Shiotani (1953) and Matsuzawa
and Takeuchi (2002) as z − 0.30
κu∗
u0 (z)min 30, 30 0.15
(3)/β.
vf (z) = α
ρpIt was noted that α and β were the constants to ensure
consistency with the period of the drifting snow eventand to relate to snowdrift density, respectively. We set
β to 0.5 for compacting snow particles, which was consistent with the condition in the site because continuous
snowfall was recorded at the nearest AMeDAS site on the
day before the drifting snow event. It had enough time toform the compacting snow layer. The friction velocity u∗
was estimated by the method in the prior paragraph. The
reference height and snow concentration in Eq. (3) referenced to Matsuzawa and Takeuchi (2002).
4.2 Simulation set‑up

Fig. 3 Schematic of incident angle θin, elevation angle θr ,and azimuth angle ϕr in the rebound process. Green spheres showsnow particles

The simulation domain covered 18 m width by 15 m
depth by 5 m height on a grid spacing of 0.05 m (Fig. 4).
The solid fence in the simulation was set at the centerof the y-axis and located at x = 0 , which was 5 m fromthe western boundary. The fence was 6 m length, 0.1 m
thickness, and 1 m height. The integration time of the
CFD calculation was 3600 s after a spin-up of 30 s
because it took about 10 s to stabilize the calculation

Tanji et al. Progress in Earth and Planetary Science(2023) 10:67

Page 5 of 12

Fig. 4 Schematic of the calculation domain and the initial wind vectors in the experiments. The black panel at x = 0 is the fence(not shown), the time step was 0.5 ms, and the result
output was every 1 s.
We calculated trajectories of a finite number of snow
particles in the snow particle module with a time interval of 1 ms until all particles were deposited. The snow
particle calculations were conducted with several initial times for the wind, and the snowdrift result was thesum of the members’ results (Fig. 5). The initial windprofile was given to the snow particle module every
second (total of 3600 members), and the wind profile
was renewed every second in the calculation. The windprofile data at 3600 s were constantly imposed when
the calculation extended beyond 3600 s. When a snow
particle fell on the surface, the snow volume defined in
Eq. (3) was released and made snowdrift on the grid. It
was noted that snow particles did not collapse and were
not eliminated by melting or sublimation processes inthis simulation. This model excluded the splash thatsnow particles pop up from the snow surface due toimpact of other snow particles collided the surface,because we only calculated limited representative particles mentioned in Sect. 4.1. According to Okaze et al.
(2018), it showed the rebound particles was muchgreater than splash particles and the rebound process
was more dominant than the splash around the surface.

5 Results
5.1 Inflow condition

Figure 6 shows the observed time series of the windspeed and snow volume flux on February 18, 2019. We
calculated the snow volume flux by the SPC outputs
which measured numbers of snow particles passing thesensor phase and their diameter. The wind speed andsnow volume flux were large between 0830 and 1300
JST (JST = UTC + 9), with a peak at about 1010 and
1040 JST, respectively. The highest wind speed was over
11 ms−1 at a height of 7 m and the largest snow volume
flux was about 8 × 10−4 m3 m−2 10min−1 at a height of
1 m in the period. The wind speed dropped to under
6 ms−1 after 1500 JST.
Figure 7 shows the vertical profiles of the observed
wind speed averaged over former 1 h and the windspeed estimated with the least squares method. The estimated wind speed at each time almost corresponded tothe observed values and increased during the driftingsnow event from 1000 to 1200 JST (Fig. 7b–d), although
there was a slight difference between the estimated andobserved wind speed at a height of 3 m. The snow volume flux was calculated by Eq. (3) with these logarithmic
wind profiles. The estimated amount of the snow volume flux overestimated the observations below a heightTanji et al. Progress in Earth and Planetary Science(2023) 10:67

Page 6 of 12

Fig. 5 Schematic of the wind data set calculated in the computational fluid dynamics (CFD) module and its treatment in the snow particle module.
Blue vectors show that we calculate the snow particles’ movements for a maximum of 70 s

Fig. 6 Time series of the observed wind (dotted lines) and snow volume flux (solid lines) at heights of (black) 1 m, (blue) 1.5 m, (pink) 3 m,and (green) 7 m from 08:00 to 16:00 JST on 18 February 2019of 4 m, especially during the drifting snow event (Fig. 8).
However, the estimated vertical profiles of the snow volume flux clearly indicated the observed feature of larger
snow volume flux at lower levels during the drifting snow

event. Therefore, we decided that the 4 h drifting snow
event occurred from 900 to 1300 JST, and α of Eq. (3) was
set to 4 because the observed snow volume flux increasedin this period (Fig. 8b–e).

Tanji et al. Progress in Earth and Planetary Science(2023) 10:67

Page 7 of 12

Fig. 7 Averaged vertical wind profiles every hour from (a) 09:00 to (h) 16:00 JST. Yellow lines show the log profiles estimated with the least squares
method, and open circles show the observed values

Fig. 8 Averaged snow volume flux profiles every hour from (a) 09:00 to (h) 16:00 JST. Yellow lines show the values estimated using Eq. (2)and the wind profiles in Fig. 7. Open circles show the observed values

Figure 9 shows the averaged vertical profiles of thewind speed and snow volume flux over the 4 h driftingsnow event in the observation. The friction velocity wasestimated as 0.297ms−1 with the least squares method.

The estimated wind profile (Fig. 9a) almost corresponded
to the observed values. The snow volume flux estimatedby the friction velocity was still larger than the observed
data at lower heights. The estimated snow volume flux

Tanji et al. Progress in Earth and Planetary Science(2023) 10:67

Page 8 of 12

Fig. 9 Vertical profiles used in the model of the wind speed and the snow volume flux averaged over the drifting snow period. Yellow lines show
the estimated profiles, and circles show the observed values

1.9 × 10−6 ­m3 ­m−2 ­s−1 was equivalent to the volume of
1.5 × 106 ­m−2 ­s−1 snow particles with 135 µm diameter at a height of 1 m. However, we considered that theestimated profile was in the range of uncertainty in theobserved data because the maximum of the observedsnow volume flux was one or two orders larger per second than the averaged observation data, especially ata height of 1 m. The estimated wind profiles were given
as the western boundary condition in the CFD module,
and the amount of the inflow snow was given in the snow
particle module.
5.2 Model results

Figure 10 shows a snapshot of the wind profiles at the endof the calculation in the CFD module. The vertical windprofile had three obvious features, which were a stagnation area formed around the front of the fence, stronger
wind speed at the top of the fence, and reverse flow
developed behind the fence from x = 0–5 m (Fig. 10a).
The reverse flow behind the fence was clearer on the horizontal wind profile (Fig. 10b). The wind flow horizontallyforked around the fence, and there was a dipole patternfrom x = 0–4 m along y = 5 and 10 m. The qualitatively
features of the vertical and horizontal wind profiles wereconsistent with previous simulation studies (Han et al.
2021; Tanji et al. 2021).
Figure 11a shows the model result of the snowdrift
distribution around the fence. The higher snowdrift wascloser to the fence on the windward side, but the heightof the snowdrift just in front of the fence was low. The
snowdrift height estimated by the model agreed with theobserved profile in the cross section along the center ofthe fence (Fig. 11b), with a peak around x = −1.5 m anda peak height of about 0.41 m. Many previous simulation
studies (Uematsu et al. 1991; Alhajraf 2004; Tanji et al.
2021) and the previous observation (Tabler 1994) also
showed this outline form of the snowdrift. This snowdrift formation was caused by the slower wind speed
near the fence (Fig. 10a). The results for the snowdriftwith half the amount of inflow particles are shown in
Fig. 11b to demonstrate the growth process of the snowdrift. The result showed that the snowdrift around the
inflow boundary mainly grew. The snowdrift nearer thefence was formed later by the snow particles that reachedthere by the repeat resuspension and rebound processes
from the mature snowdrift surface. Snowdrifts withheights of about 0.45 m were displayed from x = −0.1
–0 m only in the model, contrary to the observations witha small snowdrift just behind the fence around x = 0.5 m(Fig. 11b). The difference could be because of a small gap
between the panel and the snow surface at the observation site. Some snow particles were blown to the leeward side of the fence from the gap and accumulated justbehind the fence. This gap could cause the differences inthe snowdrift height between the simulation and observations around x = 0 m.
The snowdrift distribution on the leeward side of thefence had an arc formation in the no-fence zone, but thesnowdrift was not formed just behind the fence. Thesefeatures were consistent with previous simulation studies(Beyers et al. 2004; Okaze et al. 2013; Tanji et al. 2021).
The arc formations were made by snow particles following the split flow from the sides of the fence (Fig. 10b).
Few snow particles in the flow reached just behind theTanji et al. Progress in Earth and Planetary Science(2023) 10:67

Page 9 of 12

Fig. 10 Snapshots of the wind vectors around the fence at the last segment in the cross section along a y = 7.5 m and b z = 0.5 m. The solid lines
at x = 0 show the fences

Fig. 11 a Snowdrift height distribution around the fence. b Snowdrift height in the cross section along y = 7.5 m at (solid) 1 h and (dotted) 0.5 h
lead times. The blue circles show the observed heights

Tanji et al. Progress in Earth and Planetary Science(2023) 10:67

fence, and almost no snowdrift was formed along thecenter of the fence on the leeward side in both the simulation and observations (Fig. 11b). The peak of the snowdrift on the leeward side was y = 12.5 and 2.5 m fromx = 3 to 4 m.

6 Conclusions and discussion
In this work, we upgraded the LBM-based snowdrift
model developed by Tanji et al. (2021) with parallelization and a quantitative treatment for snow particles,
which enabled us to perform high-speed, long-run, andpractical simulation of the snowdrift. We conducted full
three-dimensional CFD and snow particle computation that reproduced snowdrifts along the center of thefence observed after one case of the drifting snow eventin the flat snowy environment at Teshikaga, Japan. In theexperiment, we used the observed vertical wind profileat the Teshikaga site as the inflow boundary conditionsin the model. The numerical simulation was conductedin a finite channel containing a finite solid fence, mimicking the observation site. The simulated result of thesnowdrift height corresponded to the observed snowdrift
height along the center of the fence. We demonstrated
that the upgraded model reproduced the observed snowdrift distribution along the center of the snow fence
quantitatively (Additional file 1).
We did not directly use the observed amount of drifting snow particles as the boundary condition due to theinsufficient vertical resolution above the snow surface.
Instead, we used the estimated value of the observed
wind data. Tanji et al. (2021) showed that the inflow
snow particles from under 1 m height above the surface are important in forming the snowdrift around the
1 m-high fence. However, measuring snow volume flux
just above the surface is difficult because the observation instruments can be buried in the snowdrift during
the events. Even the trench method, in which the number of snow particles is measured in trenches on the surface (Nishimura 2009), is not practical currently due toits high costs. Hence, the estimation of the inflow snow
particles in this paper is practical if the assumption of thelogarithmic wind profile [Eq. (2)] is valid.
The model validation in this study was done by comparison with the one-dimensional observation data forsnow accumulation along the center of the fence. Thisgauge observation of the cross section at the fence centerdid not show whether the fork structure on the leeward
side of the fence was well simulated in the model. Photogrammetry and LiDAR (light detection and ranging)have been used to measure the horizontal distribution ofsnowdrift heights. Photogrammetry mainly uses photographs taken by UAVs (unmanned aerial vehicles). Niiya
et al. (2021) took aerial images from a UAV before andPage 10 of 12

after a drifting snow event and made a digital surface
model of a snowdrift around a snow fence. The estimatedsnowdrift height was consistent with the cross-sectional
observation. Although photogrammetry is a useful tool
for observing snowdrifts, UAVs cannot fly under strong
wind conditions and images cannot be taken at night.
LiDAR can mitigate these disadvantages of photogrammetry because the method measures reflected near-infrared radiation or ultraviolet radiation. LiDAR can also
observe the snowdrift growing during a drifting snow
event. Okaze et al. (2022) showed that the snowdrift
distribution observed by LiDAR corresponded to thatobserved by photogrammetry. Ohara et al. (2022) compared a snowdrift model simulation with these observations. We will validate our model with these methods infuture work.
The numerical simulation model in this study ignored
the splash process and the sublimation process of snow
particles, which were often implemented in the modelof drifting snow in the previous studies. The splash process was considered not to contribute to the snowdrift
distribution in this situation because of more rarely
development than the rebound process, as mentioned
in Sect. 4. However, previous studies suggested that thesplash process had an important role to the development of the drifting snow (Niiya and Nishimura. 2022;
Zwaaftink et al. 2014). Therefore, the splash process isrequired when we simulate the beginning of the drifting snow event explicitly. The sublimation process also
does not have the impact on the snowdrift distribution
because the amount of snow is lost by the sublimationin the upper area of the drifting snow development. The
upper height had little snow volume flux during the drifting snow event (Figs. 8, 9). In addition, Tanji et al. (2021)
demonstrated that few snow particles inflowing fromover 1.5 m height fell on around the fence and most ofthem were blown to out of the calculation domain. The
sublimation process can be important in the discussion
about the hydrological cycle and surface water mass balance during the entire winter season in the cryosphere(Zwaaftink et al. 2013; Palm et al. 2017).
The present model in this study remained some room
for improvement, but the model had the potential to be
applied to resolving many problems caused by snowdrifts around roads. Snow fences have been installed on
the side of the roads to protect the roads from snowdrifts on the leeward. However, no one knows whichsnow fence structures are optimal for preventing snowdrifts developing on roads. The LBM snowdrift model
is suitable for researching structures in many experiments because the model has high parallel computation efficiency. Snowdrifts that form around the lateralgap between snow fences are also a serious problem forTanji et al. Progress in Earth and Planetary Science(2023) 10:67

road administration because gaps are required at intersections. Snowdrifts around the various road structures
can be estimated with this model. This model also canbe used to estimate visibility due to drifting snow onroads, which is the main cause of traffic disruption insnowy areas. Visibility due to drifting snow is related
to the amount of snowdrift flux (Budd et al. 1966;
Takeuchi and Fukuzawa 1976; Matsuzawa and Takeuchi2002). The amount of drift snow transport can be calculated explicitly in this model by tracing drifting snow
particles.
In mountainous areas, the redistribution of snow particles by drifting snow plays an important role in snow
cover distribution and contributes to the avalanche
potential, but previous studies did not calculate theredistribution explicitly (Bartelt and Lehning 2002; Brunet al. 1992). The boundary conditions and snowdrift distribution in mountainous areas are more complicatedthan those around a snow fence, although the LBM canstill calculate the turbulent flow because it is suitablefor complicated boundary conditions. Tracing driftingsnow particles is also helpful for estimating the avalanche
potential because the potential is strongly associated withthe condition of the snowpack layers. The snowpack layer
consisting of drifting snow particles has different featuresfrom that consisting of snowfall particles because theshape of the particles is different. Therefore, the model inthis paper has the potential to resolve problems causedby drifting snow in mountainous areas.
Abbreviations
AMeDAS	Automated Meteorological Data Acquisition System
CFD	Computational fluid dynamics
LBM	Lattice Boltzmann method
LES	Large-eddy simulation
LiDAR	Light detection and ranging
SPC	Snow particle counter
UAV	Unmanned aerial vehicle

Supplementary Information
The online version contains supplementary material available at https://​doi.​
org/​10.​1186/​s40645-​023-​00599-3.
Additional file 1. Supplementary Material.
Acknowledgements
We thank Dr. Yousuke Sato for teaching basic parallelization. We are very grateful to Dr. Satoshi Omiya for supporting us in the observations. We also thankDr. Tsubasa Okaze, Dr. Hirofumi Niiya, Dr. Masaki Nemoto, and Prof. Kouichi
Nishimura for providing insightful comments in the early stage of this study.
Author contributions
ST, the main author, contributed to this paper, conceived and designed thestudy, developed the source code of the snowdrift model, analyzed the model
results, and wrote the manuscript. MI, the second author, helped in the development of the snowdrift model, interpretation of the results, and writing themanuscript. YH contributed to observation work and acquired and analyzed
the observed data in this study.

Page 11 of 12

Funding
ST was supported by JST SPRING, Grant Number JPMJSP2119. MI is supported
by JSPS KAKENHI Grant No. 19H00963, by the Environment Research andTechnology Development Fund JPMEERF20232003 of the Environmental Restoration and Conservation Agency provided by Ministry of the Environmentof Japan, and by Research Field of Hokkaido Weather Forecast and Technology Development (endowed by Hokkaido Weather Technology Center Co.,
Ltd.). MI and YH were supported by grant number JPMXD0722680734 from
the Ministry of Education, Sports, Culture, Science and Technology (MEXT).
This work was partly performed by using the supercomputer system at theInformation Initiative Center, Hokkaido University, Sapporo, Japan.
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.
Received: 2 June 2023 Accepted: 9 November 2023

References
Alhajraf S (2004) Computational fluid dynamic modeling of drifting particles atporous fences. Env Modelling Softw 19:163–170. https://​doi.​org/​10.​1016/​
S1364-​8152(03)​00118-X
Bartelt P, Lehning M (2002) A physical SNOWPACK model for the Swiss avalanche warning. Part I. numerical model. Cold Reg Sci Tech 35(3):123–145.
https://​doi.​org/​10.​1016/​s0165-​232x(02)​00074-5
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
Brun E, David P, Sudul M, Brunot G (1992) A numerical-model to simulate
snow-cover stratigraphy for operational avalanche forecasting. J Glaciol
38(128):13–22. https://​doi.​org/​10.​3198/​1992J​oG38-​128-​13-​22
Budd WF, Dingle R, Radok U (1966) The byrd snow drift project: outline andbasic results. Studies in antarctic meteorology. Antarct Res Ser 9:71–134.
https://​doi.​org/​10.​1029/​AR009​p0071
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
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
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 Meth Engng 72:1111–1134. https://​doi.​org/​10.​1002/​nme.​2114
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
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 P01018.
https://​doi.​org/​10.​1088/​1742-​5468/​2010/​01/​P01018.
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 Meteorol
164:1651–2181. https://​doi.​org/​10.​1007/​s10546-​017-​0249-y
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

Tanji et al. Progress in Earth and Planetary Science(2023) 10:67

Matsuzawa M, Takeuchi M (2002) A study of methods to estimate visibility
based on weather conditions. Seppyo 64:77–85. https://​doi.​org/​10.​5331/​
SEPPYO.​64.​77 (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
Niiya H, Nishimura K (2022) Hysteresis and surface shear stresses during snowparticle aeolian transportation. Bound-Layer Meteorol 183:447–467.
https://​doi.​org/​10.​1007/​s10546-​022-​00688-8
Niiya H, Omiya S, Sunako S, Nishimura K, Okaze T (2021) Snowdrift observation
around fence by UAV-SfM photogrammetry. JSSI JSSE Joint Conference
2021 in Online (in Japanese)
Nishimura K (2009) Measurement of Blowing Snow. Nagare 28:455–460 (inJapanese)
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
Nishio F, Ishida T (1973) Rate of Turbulent Energy Dissipation During Snow
Drifting. Low Temp. Sci. Series A Phys Sci 31:69–85 (in Japanese)
Ohara N, He S, Parsekian AD, Jones BM, Rangel RC, Nichols IO, Hinkel KM (2022)
Spatial snowdrift modelling for an open natural terrain using a physicallybased linear particle distribution equation. Hydrol Processes 36:1. https://​doi.​org/​10.​1002/​hyp.​14468
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
Okaze T, Niiya H, Omiya S, Sunako S, Nishimura K (2019) Observation of snowdrift around two-dimensional fence during a drifting snow event. JSSI &JSSE Joint Conference 2019 in Yamagata 110 (in Japanese)
Okaze T, Niiya H, Omiya S, Sunako S, Kawashima R, Tanji S, Nishimura K (2022)
Measurement of Snowdrift around a Fence with LiDAR for In-vehicle Use.
JSSI & JSSE Joint Conference 2022 in Sapporo (in Japanese)
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. https://​doi.​org/​10.​3130/​aije.​78.​149 (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 area inTokyo. High Perf. Comput. Symp. 123–131 (in Japanese)
Palm SP, Kayetha V, Yang Y, Pauly R (2017) Blowing snow sublimation and transport over Antarctica from 11 years of CALIPSO observations. Cryosphere11(6):2555–2569. https://​doi.​org/​10.​5194/​tc-​11-​2555-​2017
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 Wea 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. Comp & Math
App 69:518–529. https://​doi.​org/​10.​1016/j.​camwa.​2015.​01.​010
Sugiura K, Maeno N (2000) Wind-tunnel measurements of restitution coefficients and ejection number of snow particles in drifting snow: Determination of splash functions. Bound-Layer Meteorol 95:123–143
Tabler RD (1994) Design guidelines for the control of blowing and driftingsnow. Strategic Highway Research Program Report SHRP-H-381.
Takeuchi M, Fukuzawa Y (1976) On the light attenuation and visibility in snow
drift. Seppyo 38(4):165–170. https://​doi.​org/​10.​5331/​seppyo.​38.4_​165 (inJapanese)Tanji S, Inatsu M, Okaze T (2021) Development of a snowdrift model with thelattice Boltzmann method. Prog Earth Planet Sci 8:57. https://​doi.​org/​10.​
1186/​s40645-​021-​00449-0
Uematsu T, Nakata T, Takeuchi K, Arisawa Y, Kaneda Y (1991) Three-dimensional
numerical simulation of snowdrift. Cold Reg Sci Tech 20:65–73. https://​doi.​org/​10.​1016/​0165-​232X(91)​90057-N
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 Problem
Eng 742432. doi:https://​doi.​org/​10.​1155/​2014/​742432.
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

Page 12 of 12

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, US, p 640
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
Zwaaftink CDG, Mott R, Lehning M (2013) Seasonal simulation of drifting snow
sublimation in Alpine terrain. Water Resour Res 49:1581–1590. https://​doi.​
org/​10.​1002/​wrcr.​20137

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