The Visual Computer manuscript No.
(will be inserted by the editor)

Visual Simulation of Mixed-motion Avalanches with Interactions
Between Snow Layers
Yusuke Tsuda · Yonghao Yue · Yoshinori Dobashi · Tomoyuki Nishita

Abstract In the ﬁeld of computer graphics, simulation of
ﬂuids, including avalanches, is an important research topic.
In this paper, we propose a method to simulate a kind ofavalanche, mixed-motion avalanche, which is usually large
and travels down the slope in fast speed, often resulting in
impressive visual effects. The mixed-motion avala- nche consists of snow smokes and liqueﬁed snow which form an upper suspension layer and a lower dense-ﬂow layer, respectively. The mixed-motion avalanche travels down the surface of the snow-covered mountain, which is called accumulated snow layer. We simulate a mixed-motion avalanchetaking into account these three snow layers. We simulate
the suspension layer using a grid-based approach, the denseﬂow and accumulated snow layer using a particle-based approach. An important contribution of our method is an interaction model between these snow layers that enables us toobtain the characteristic motions of avalanches, such as thegeneration of the snow smoke from the head of the avalanche.
Keywords Simulation · Fluids · Avalanche · Generation ofsnow smoke · Adhesion of snow smoke · Snow entrainment

1 Introduction
In the ﬁeld of computer graphics, simulation of natural phenomena is one of the most important research topics. Among
various natural phenomena, avalanches need to be simulated
by computer graphics because of the difﬁculties to take photographed images and videos.
Y. Tsuda · Y. Yue · T. Nishita
The University of Tokyo
E-mail: {tsuday, yonghao, nis}@nis-lab.is.s.u-tokyo.ac.jp
Y. Dobashi
Hokkaido University
E-mail: doba@ime.ist.hokudai.ac.jp

Flowing avalanche

Powder-snow avalanche

Mixed-motion avalanche
Snow smoke

Liquefied snow(a)Snow smoke(b)Liquefied snow(c)Fig. 1 (a) ﬂowing avalanche, (b) powder-snow avalanche, and (c)mixed-motion avalanche.

Avalanches are classiﬁed into three categories: 1) ﬂowing avalanches, 2) powder-snow avalanches and 3) mixedmotion avalanches [10] (see Figure 1). In this paper, we propose a method to simulate the mixed-motion avalanche.
The mixed-motion avalanche consists of two-tiered snow
layers, called dense-ﬂow layer and suspension layer, shownas Figure 1 (c). The dense-ﬂow layer is the lower layer consisting of the snow which behaves like a liquid (such snow
is called liqueﬁed snow). The suspension layer is the upper layer consisting of snow smoke and air, covering thedense-ﬂow layer. The snow on the surface of a snow-covered
mountain, called accumulated snow layer, is also important for the motion of an avalanche, because the dense-ﬂow
layer usually entrains an accumulated snow layer, making
the avalanche grow larger during the motion down the slope.
Interactions between these three snow layers are important
to obtain the motion of avalanche. An important contributionof this paper is an interaction model for the mixed motion
avalanche between these snow layers.
In our method, the snow in the suspension and denseﬂow layers are simulated as ﬂuids. The dense-ﬂow layer isabout 100 times denser than the suspension layer [10]. We
simulate the suspension and dense-ﬂow layers as two-phaseincompressible ﬂuids for the stability of the simulation. The
snow of the suspension layer and that of the dense-ﬂow layer
can be assumed to be the ﬂuids with relatively small density
differences, respectively [10]. This assumption is supportedby the observations of real avalanches, as described in theliterature [10].

2

The snow smoke in the suspension layer is simulated
by a grid-based approach [16], because a number of particles are necessary to calculate the detail of the snow smoke[4]. On the other hand, the dense-ﬂow layer is simulatedby Smoothed Particle Hydrodynamics, SPH [5][12], because
the dense-ﬂow layer may undergo large deformations, suchas the splash of the snow of the dense-ﬂow layer when thedense-ﬂow layer runs into obstacles.
In our method, we take into account two types of interactions: 1) interactions between the suspension and denseﬂow layers and 2) interactions between the dense-ﬂow andaccumulated snow layers. Firstly, we assume that there isan overlapping region with a high density contrast between
the dense-ﬂow and suspension layers [10]. We assume thatthe overlapping region has non null thickness, and both thesnow smoke and liqueﬁed snow exist and interact with eachother in the overlapping region. Secondly, we take into account snow entrainment: the dense-ﬂow layer entrains thesnow of the accumulated snow layer.
We show that our method can reproduce the characteristic phenomena of the mixed-motion avalanches observedin the real world, such as the generation of the snow smokefrom the head of the avalanche.
2 Related Workels of recent years and depth-averaged models, and those of
the powder-snow avalanches can be calculated by density
current models. However, all of these four models cannot
be used to simulate the mixed-motion avalanches because
both the suspension and dense-ﬂow layers are not taken into
account. In our method, both these layers are taken into account, together with an interaction model between snow layers. The dense-ﬂow layer is calculated by using particles asin the granular ﬂow models. The suspension layer is calculated by using grid taking into account the inﬂuence of the
air ﬂow to the motion of the snow smoke, as in the density
current models.
3 Previous Methods for Fluid Simulation
In this section, we brieﬂy review grid-based and particlebased ﬂuid simulation methods. The motions of the ﬂuids
are simulated by solving incompressible Navier-Stokes Equations:∇ · u = 0,(∂uρ+ (u · ∇)u
∂t

(1)
)
= −∇p + µ ∇2 u + f,(2)

where u is the velocity, t is time, p is the pressure, ρ is theﬂuid density, µ is the viscosity and f is external forces.

In the ﬁeld of computer graphics, simulation models for theavalanches based on the physical properties of the snow lay3.1 Grid-based Method
ers have not been proposed to our knowledge. Kapler [8]created a computer-generated animation of an avalanche covIn grid-based method [16], the quantities of the ﬂuids, such
ered by snow smoke. In [8], snow smoke is represented byas the velocity and density, are stored at ﬁxed grid points,
3D volume texture attached to a series of rotating particles
and the inﬂows and outﬂows of the quantities are calculated
emitted from the head of an avalanche. In their results, theat each grid point. The velocities on grid points are updated
head of the avalanche moves through the air and emits thesnow smoke in the air even though the head of the real avalanche by calculating Eq. (2), and a scalar quantitiy a is updated
according to:usually moves on the ground and emit the snow smoke into
the air.
∂a
In other research ﬁelds, i.e., glaciology, physics, and com= −(u · ∇)a + κ ∇2 a + S,
(3)∂tputational geoscience, simulation models of avalanches havewhere κ is the diffuse rate of the quantity and S is the inbeen proposed. However, no universal dynamic model has
crease or decrease of the quantity due to the external source.
been established, and each proposed dynamic model is basedon assumptions depending on categories of the avalanches.
Previous simulation models for avalanches can be classiﬁed
3.2 Particle-based Method
into following four types: center-of-mass models [13], granular ﬂow models [14], depth-averaged models [7] and denIn SPH [12], the motions of ﬂuids are represented by thesity current models [2].
movement of particles. Each particle stores the position, velocity and mass. The motion of the ﬂuid is simulated by folAmong these four models, center-of-mass models andlowing steps: 1) calculation of the densities of the ﬂuid at
old granular ﬂow models cannot be used to calculate thethe center of the particles, 2) calculation of the forces acting
shape and deformation of the avalanche, because the centeron particles, 3) update of the velocities and positions of theof-mass models and old granular ﬂow models regard an avalparticles according to the calculated forces.
anche as a simple mass point and a continuous body without
The density, pressure force and viscosity force of a pardeformation, respectively. The shape and deformation of theticle are calculated based on the weighted inﬂuences on theﬂowing avalanches can be calculated by granular ﬂow mod-

3
Snow layers
Suspension layer(Grid)
Dense-flow layer(Particles)
Accumulated snow layer(Particles)

Interactions between
snow layers

Simulation of dense-flow layer

Between suspension layerand dense-flow layer
Between dense-flow layer
and accumulated snow layer

Fig. 2 Simulation model of a mixed-motion avalanche in the ourmethod: Green mesh indicates a grid used in grid-based method. Blueand red circles indicate the SPH particles for the dense-ﬂow layer andthe accumulated snow layer, respectively. Orange arrows indicate interactions between dense-ﬂow layer and suspension layer. The light green
region indicates the overlapping region of the suspension and denseﬂow layers. A part of the particles in the accumulated snow layer has
moved from default position due to snow entrainment.

particle from neighboring particles. In SPH, a scalar or vector quantity A(r) at the position r in the ﬂuid is calculated
as the weighted sum of the quantity Ai of the neighboring
particle i:()Ai
A(r) = ∑ mi W |r − ri |, h ,ρiiInteractions between suspension
and dense-flow layers

Previoustime step

Simulation of suspension layer

Interaction between dense-flow
and accumulated snow layers

Nexttime step

Fig. 3 Overview of our method: Processes in the ﬁgure are calculated
in one time step. 1) Interaction forces between particles and grid points
are calculated. 2) The motion of the dense-ﬂow layer is calculated by
SPH. 3) The motion of the suspension layer is calculated by a gridbased approach. 4) Snow entrainment is calculated.

(4)

where mi is the mass of the particle i, ρi is the density of theparticle i, W is a smoothing kernel function calculating theweight of inﬂuence on a particle from a neighboring particle, according to the distance between these two particles. h
is used as the radius of support for weight calculation, andits value is constant and is same for all the particles. The
density, pressure force and viscosity force of particle i arecalculated by using Eq. (4). The smoothing kernels used to
calculate these quantities are introduced in [12].
4 Proposed Method
Figure 2 shows the simulation model of our method. We
simulate the motion of the suspension layer using a gridbased approach, and the dense-ﬂow and the accumulated
snow layers using a particle-based approach. Additionally,
we introduce an interaction model between the snow layers.
We take into account the interactions between the suspension layer and the dense-ﬂow layer, and the interactions between the dense-ﬂow layer and the accumulated snow layer.
We ignore the interaction between suspension layer and theaccumulated snow layer. The interactions between the suspension layer and the dense-ﬂow layer are assumed to occur
in the overlapping region (described in detail in Section 4.1),
and such interactions include the generation and adhesion of
the snow smoke and the interaction forces acting in the overlapping region. As the interactions between the dense-ﬂow
layer and the accumulated snow layer, we consider the snow
entrainment. Figure 3 shows the ﬂowchart of the processes
calculated in a single time step.
In the following sections, we ﬁrst deﬁne the overlapping
region in Section 4.1. Then, we describe the details of thesimulation of the suspension layer in Section 4.2, and denseﬂow and accumulated snow layers in Section 4.3. We describe the interactions between snow layers in Section 4.4.
Finally, we describe the initial conditions to start the simulation in Section 4.5.

4.1 Deﬁnition of the Overlapping Region
Since there is no contact surface between the suspension anddense-ﬂow layers, we do not use previous methods for themulti-phase ﬂuids, such as [6], which are designed to simulate the shapes and motions of the contact surfaces of the
ﬂuids. Instead, we introduce a concept of overlapping regions, in which both of the snow smoke of the suspension
layer and the liqueﬁed snow of the dense-ﬂow layer exist,
and the snow of these two layers interact with each other.
The overlapping regions are deﬁned as the region where thenumber densities of the particles are relatively low comparedto those in the dense-ﬂow layer. The number density n(x) at
position x is calculated by Eqs. (5) and (6) [9].
()n(x) = ∑ Wn |r j − x|, h ,
(5)
j

{h
Wn (r, h) =r − 1 (0≤r < h)0(r > h),(6)where Wn is a kernel function for particle number density, r
is the distance between particles j and position x, and h is theradius of support for weight calculation. If n(x) is positiveand lower than a threshold n0 , we regard position x to be in
the overlapping region. Figure 4 shows an example of theclassiﬁcation of the particles in the dense ﬂow layer.

4Fig. 4 Classiﬁcation of the particles. The particles of the dense-ﬂow
layer traveling down the slope are shown without the snow smoke ofthe suspension layer. Particles in the overlapping region are shown inpurple color, and the other particles are shown in yellow color.

4.2 Details of the Simulation of Suspension Layer
In our method, each grid point stores the velocity, the ﬂuid
density and the ratio rs of the snow smoke density to thetotal density:rs =ρsmoke,ρair + ρsmoke(7)

where ρair is the density of the air, and ρsmoke is the densityof the snow smoke. The ﬂuid density ρ in Eq. (2) is the total
density of the snow smoke and the air, i.e., ρ = ρsmoke + ρair .
We substitute rs for a in Eq. (3) when solving Eq. (3). External forces acting on grid points are calculated by Eq. (8).
f = ρ g + fbuoyancy + fdrag + fli f t + fcon f ,(8)where g is a gravitational vector, fdrag and fli f t indicate theinteraction forces on a grid point from neighboring particles. fdrag is the force that makes particles move along with
the ﬂow on the neighboring grid points, and acts as an air
resistance force on the particle. fli f t is the force that makes
the particle move according to the vortices (see Section 4.4.1
for details). fbuoyancy is the buoyancy acting on snow smoke,and is calculated by Eq. (9).
fbuoyancy = −(1 − rs )ρ g.

Fig. 5 Yellow line indicates the boundary for the suspension layer, andlight green region indicates the overlapping region. Purple points indicate the grid points determined to be out of the boundary, and grid
points on the yellow line are the grid points on the boundary.

Secondly, we treat a grid point as on the boundary when thegrid point is not out of the boundary and one or more adjacent grid points are out of the boundary. We let the boundary
condition be that the normal component of the velocity to theboundary is zero.

4.3 Details of the Simulations of Dense-Flow and
Accumulated Snow Layers
In our method, each SPH particle stores its position, velocity and mass of the snow. As external forces, a gravitational
force and reaction forces of fdrag and fli f t in Eq. (8) from theneighboring grid points are applied to the particles.
The particles of the accumulated snow layer are not allowed to move initially, and begin to move when the particles meet the condition for snow entrainment (see Section
4.4.2 for details). After a particle of the accumulated snow
layer begins to move, the particle is treated as a particle ofthe dense-ﬂow layer.
Figure 3 shows the ﬂowchart of the processes calculated
in a single time step.
4.4 Interactions between Snow Layers(9)fcon f is an external force used to generate vortices calculated
by vorticity conﬁnement [3].
Next, we describe the boundary condition for the suspension layer. The boundary for the suspension layer is deﬁned as the boundary between the accumulated snow layerand suspension layer, and between the overlapping region
and the region of the dense-ﬂow layer out of the overlapping region (see Figure 5). We ﬁnd grid points that are on
the boundary, as follows. Firstly, we ﬁnd grid points that are
obviously out of the boundary. These grid points are the ones
that meet at least one of the following conditions:(1) The grid point is the nearest to more than one particlesof the dense-ﬂow layer, and all of these particles are in
the overlapping region.
(2) The grid point is not the nearest grid point to any of theparticles.

Firstly, we introduce interactions between the dense-ﬂow
layer and suspension layer: generation of the snow smoke,adhesion of the snow smoke and interaction forces. In theoverlapping region, the snow of the dense-ﬂow layer is blown
up and becomes snow smoke. On the other hand, the snow
smoke adheres to the dense-ﬂow layer. In this paper, we callthese phenomena “generation of the snow smoke” and “adhesion of the snow smoke”, respectively. Addition to theseinteractions, the ﬂuids in the suspension and dense-ﬂow layers push each other in the overlapping region through theinteraction forces.
Secondly, we introduce snow entrainment as the interaction between the dense-ﬂow and accumulated snow layers.
The snow of the accumulated snow layer is entrained by thatof the dense-ﬂow layer and the total mass of the avalanche
increases. Snow entrainment inﬂuences the size and speedof the avalanches and the deformation of the surface of theaccumulated snow layer.

5Fig. 6 The snow smoke is generated from the particle according to therelative velocity, indicated by the red arrow. The snow smoke is shownby blending white color with the background.

4.4.1 Interactions between Suspension and Dense-ﬂow
LayersAs the forces between particle and grid point, we used coupling forces as described in Eqs. (10) and (11) [11].
fdrag
i, j = −kdragmi|vi − ui |(vi − ui ),riflii, jf t = −kli f t mi (vi − ui ) × ωi ,to the particle i, using Eq. (13), based on the intuition thatthe amount of adhesion is related to the swept volume.
Ai = ∑ Ca d j Sisec |u j − vi |∆ t,

(13)

j

(10)
(11)

li f twhere fdragi, j and fi, j are the forces acting on particle i fromgrid point j, and kdrag and kli f t are constants. As the values of kdrag and kli f t , we used 0.81 and 0.2, respectively.
mi , ri and vi are the mass, radius and velocity of the particle
i. ui is the velocity at the position of the particle i interpolated from neighboring grid points. ωi is the vorticity at theparticle i. ωi is interpolated from the vorticities of neighboring grid points, and vorticity of grid point j is calculated byω gj = ∇ × ugj , where ugj is the velocity of the grid point j.
Forces acting on a grid point from particles are calculated asli f tthe reaction forces of fdragi, j and fi, j .
To simulate generation of the snow smoke, a part of massof the particle in the overlapping region is transferred to
the neighboring grid point. Although the size of a particle
should change according to the mass transfer, we currently
ﬁx the size of the particle and change its mass, instead. The
amount of the transferred snow is related to the relative velocity of the particle to the grid point as shown in Figure
6. According to a real-world measurement [17] of the relationship between the wind velocity and the density blown
from the surface of the accumulated snow layer, we model
the amount G j of generated snow smoke on the grid point j
Eq. (12).

G j = ∑ Cg |u j − vi |1.7 Sisur ∆ t,

Fig. 7 The movement of the particle and adhesion progress is shown
from left to right. Red arrow indicates the movement of the particle
according to the relative velocity of the particle to the grid point. Diagonal lines indicate the sectional area of the particle.

(12)

where Ca is a user-speciﬁed constant, and we used 0.2. d j isthe density of the snow smoke on the grid point j neighborto the particle i, and Sisec is the sectional area of the particle
i. In Eq. (15), Sisec |u j − vi |∆ t indicates the volume swept by
the particle in a single time step.
4.4.2 Interactions between Dense-ﬂow and Accumulated
Snow Layers
As the interactions between dense-ﬂow and accumulated snow
layers, we take into account snow entrainment. We ignore
the effect that snow in the dense-ﬂow sticks onto the accumulated snow layer. In the snow entrainment process, theparticles of the dense-ﬂow layer entrain the particles fromthe accumulated snow layer. We determine whether eachparticle of the accumulated snow layer is entrained or not,according to the result of the experiment on the fracture of
the accumulated snow [18]. In [18], relationship between
fracture of the accumulated snow and the strain rate (deﬁned as the change in strain over the time) of the snow isresearched, and threshold of strain rate for the fracture of
the accumulated snow is obtained.
In our method, particles of the accumulated snow layer
are entrained when the calculated strain rate of each particle is greater than a threshold γ0 obtained from [18], whichranges from 0.000001 to 0.0001 exponentially to the depth
at the particle from the surface of the accumulated snow
layer. The strain rate γi is calculated based on the relationship between strain rate of snow and the mass the snow issustaining [10]:iwhere Cg is a constant, u j is the velocity on the grid point j，vi is the velocity of the particle i neighbor to the grid pointj，and Sisur is the surface area of the particle i. The user can
specify the amount of the snow smoke by determining thevalue of Cg . In our method, we used 0.2 as the value of Cg
according to [17].
We model the adhesion of the snow smoke based on theidea that each particle in the overlapping region sweeps theregion, and the snow smoke in the swept region adheres tothe particle. We calculate the amount Ai of the snow adheredγi =1
Wi sinθ + γ f ,ηη = 8.5 × 106 exp(0.021ρa ),

(14)
(15)

where η is the viscosity coefﬁcient of the accumulated snow,
and Wi is the total mass of the snow of the dense-ﬂow andaccumulated snow layers on the top of particle i [10] (seeFigure 8). We ignore the masses of the air and the snow of
the suspension layer because they are negligible. θ is the angle of the slope, and γ f is the strain rate due to friction. ρa

6

(a)
(b)
Fig. 8 Blue and red circles indicate the particles of the dense-ﬂow
layer and the accumulated snow layer, respectively. A particle of theaccumulated snow layer receives the total mass of the snow in the lightbrown region. The accumulated snow layer receives more mass of thesnow in the vertical direction when the slope is (a) more steeply inclined than (b) gently inclined.

is the density of the accumulated snow layer. As the valueof ρa , we used observational value (100kg · m−3 ). Eq. (15)
is obtained from measurements [10] and indicates the relationship between the viscosity coefﬁcient and density of theaccumulated snow. The value of η is usually large and γi isclose to the threshold γ0 .
As described above, the determination of the snow entrainment whether a particle i of the accumulated snow layer
is entrained or not is based on the comparison between γiand γ0 . However, if we conduct this determination at eachtime step, the result will be affected by the time step ∆ t:smaller ∆ t will result in more times of the determination andmore snow will be entrained when ∆ t is smaller. To resolve
this problem, we conduct the determination not at every time
step but in a certain probability Pe at each time step, whichis calculated as:
Pe =Er· ∆ t · Ssec ,
m̄

(16)

where Er is the observational mass value of the entrainedsnow from a 1m2 surface of the accumulated snow layer per
one second, 10(kg · m−2 · s−1 ) [15]. m̄ is the initial mass ofa single particle of accumulated snow layer. Ssec is the sectional area of the particle. Strictly speaking, Pe is not linear
to ∆ t, but can be approximated as Eq. (16) if Pe is sufﬁciently small, as shown in Appendix A. Each particle of theaccumulated snow layer is treated as a particle of the denseﬂow layer after the particle is entrained.

4.5 Initial Conditions
The placement of the particles in the accumulated snow layer
should be arranged in order not to make the simulation unstable when a particle of the accumulated snow layer is entrained: if the distance between two particles are too small, astrong force acts between two particles and makes the simulation unstable. We arrange the particles of the accumulated
snow layer to meet two conditions: 1) the particles are arranged nearly uniformly, 2) the average distance between
two nearest neighbor particles is the same as the value ofthe initial particle distance of the dense-ﬂow layer particles. To meet these conditions, we ﬁrst determine the initial positions of the particles by using pseudo random numbers so that these particles are arranged nearly uniformly.
Then, these particles are iteratively moved to be arranged
more uniformly by using a relaxation method [19].
We set the velocities and densities on grid points to zero(which can be set arbitrary according to user’s demand). The
particles of the dens-ﬂow layer are set on the top of the slope.
5 Results
Section 5.1 shows the results of 2D simulation, and Section
5.2 shows the results of 3D simulation. Simulations are conducted on a PC with Intel Core 2 DUO 3.33 GHz CPU. ∆ t isdetermined in an adaptive way according to CFL condition[1].
5.1 Results of 2D Simulation
In 2D simulation, a vertical cut of an avalanche is calculated.
The resolution of the grid is 690 × 240. The number of particles in the dense-ﬂow layer increased from 1,000 to 2,800
on an average, due to the snow entrainment. Memory usage
is about 40MB. Grid spacing is 0.33m, and the inﬂuential
radius of the particles is 0.21m. Each simulation consists ofabout 2,400 time steps, and the average computational time
for one time step is about 80 ms.
In the results shown in Figures 9, 10, 11 and 12, particles of the dense-ﬂow layer are shown with white color,
and the snow smoke is shown by blending white color withthe background color according to the ratio of snow smoke
density ρs on a grid point. The simulation time is shown in
the upper left corner of each result image, and the scale oflength is shown in the lower right.
Figure 9 (a) shows the result of the simulation of anavalanche traveling down a bumpy slope. The number ofthe particles of the dense-ﬂow layer increased from 1,000
to 2,560 by snow entrainment. Characteristic phenomenaof the mixed-motion avalanches that are often observed in
the real world, such as the rounded shape of the section of
the avalanche head, long tail of the avalanche, and swirlingsnow smoke are reproduced simultaneously. Figures 9 (b)and (c) show the simulation results of avalanches travelingdown the steeply and gently inclined slope, respectively. Bysnow entrainment, the number of particles of the dense-ﬂow
layer increased from 1,000 to 3,700 in (b), and from 1,000
to 2,040 in (c). More particles are entrained from the accumulated snow layer when the slope is inclined more steeply
because the accumulated snow layer receives more weightof the snow (see Figure 8). The avalanche on the steeply inclined slope travels down the slope at high speed. Therefore,
the air resistance force on the dense-ﬂow layer is strong and

7
1.0s

1.0s

1.5s

0

0

20m

3.5s

20m

3.0s

3.0s

3.0s

0

20m

3.0s

0

20m

(a)
0

20m

0

20m

0

0

20m

0

20m

(b)

20m

(a)
(b)
(c)
Fig. 9 Avalanches travelling down differently inclined slopes.
1.5s

1.5s

1.5s

0

0

0

20m

3.5s

20m

3.5s

0

20m

0

20m

3.5s

20m

(c)
(d)
Fig. 11 (a) The result of the simulation with snow entrainment. (b)The result of the simulation without snow entrainment. (c) and (d) The
accumulated snow layers after the avalanches in (a) and (b) passed.
2.5s

0

20m

0

20m

(a)
(b)
(c)
Fig. 10 Comparison of the effects of the interactions. (a) Simulation
with generation, adhesion of snow smoke and interaction forces. (b)Simulation with interaction forces but without generation and adhesionof snow smoke. (c) Simulation with generation and adhesion of snow
smoke but without interaction forces.

4.5s

0

20m

0

20m

Fig. 12 An avalanche falling off the cliff in 2D simulation.

the amount of the snow smoke generation is large, resulting
in the thick avalanches, as shown in (b). On the other hand,the avalanche on the gently inclined slope travels down theslope at low speed, resulting in the thin snow layers because of weak air resistance force and small amount of snow
smoke generation.
Figure 10 shows the comparison of the effects of theinteractions between the suspension and dense-ﬂow layers.
Figure 10 (a) is the result of the simulation with generation,adhesion of snow smoke and interaction forces. Figure 10(b) shows the result of the simulation with interaction forces
but without the generation and adhesion of the snow smoke.
The characteristic motion of the snow smoke, such as thegeneration from the head of the avalanche and swirling motions, are not seen in (b). Therefore, the simulation without taking into account snow smoke is not suitable for themixed-motion avalanche. Figure 10 (c) shows the result of
the simulation with generation and adhesion of snow smoke
but without interaction forces. Compared to Figure 10 (a),the shape of the snow smoke of the suspension layer is too
simple. Therefore, the interaction forces are also important.
In Figure 11, the results of the simulations with (a) andwithout (b) snow entrainment process are shown. In (a) and(b), the particles of the dense-ﬂow layer at the same simulation time are shown without snow smoke. White particles
indicate the particles of the dense-ﬂow layer when the sim-

ulation starts, and the dark gray particles indicate the particles entrained from the accumulated snow layer. As shownin (a), the avalanche becomes larger by snow entrainment.
The avalanche in (a) travels down the slope faster than theavalanche in (b), which agrees with the fact that larger avalanche travels down the slope faster [10]. (c) and (d) show theaccumulated snow layers after the avalanches of (a) and (b)have traveled down the slopes, respectively. The surface of
the accumulated snow layer in (c) is smoother than that in(d) because the surface snow of the accumulated snow layer
is entrained in (c).
Figure 12 shows the result of the simulation of an avalanche falling off a cliff. During falling, the snow smoke is generated from the head of the avalanche, and the head does not
spread because of air resistance. The snow in the head of theavalanche is then pushed by the snow falling from the cliff.

5.2 Results of 3D Simulation
In 3D simulation, the resolution of grid is 30 × 40 × 170,
and the number of particles of the dense-ﬂow layer is increased from 10,000 to 22,000 by snow entrainment. Grid
spacing is 1.0m, and inﬂuential radius of the particles is
0.21m. Memory usage is about 200 MB. Whole simulation

8

Fig. 13 An avalanche traveling down the slope in 3D simulation.

consists of about 3,000 time steps, and the average computational time for a single time step is about 8 seconds.
Figure 13 shows the result of our method in 3D. An
avalanche is traveling down a slope of a mountain. Characteristic phenomena of the mixed-motion avalanche, suchas the generation of the snow smoke from the head of theavalanche, are reproduced.
6 Conclusions and Future WorkIn this paper, we have presented a simulation method formixed-motion avalanches taking into account snow layers
and interactions between snow layers. The suspension layer
is simulated by a grid-based approach, and the dense-ﬂow
and accumulated snow layers are simulated by a particlebased approach. We calculated interactions between the suspension and dense-ﬂow layers and snow entrainment. We
showed that our method can simulate the characteristic phenomena of the mixed-motion avalanches observed in the realworld, such as the generation of the snow smoke from thehead of the avalanche.
As future work, we would like to use rendering techniques taking into account light scattering by snow for realistic rendering of avalanches. Moreover, we would like toreduce computational time by using GPU and more efﬁcient
representation for the suspension layer, e.g., adaptive grids.
Acknowledgements We gratefully acknowledge helpful comments andthoughtful suggestions from the anonymous reviewers. We would also
like to thank Dina Systems Co., Ltd. for providing the terrain dataof Mt. Myoko for 3D simulation. The second author is supported by
Grant-in-Aid for JSPS Fellows (20・7968).

References
1. R. Courant, K. Friedrichs, and H. Lewy. On the Partial Difference Equations of Mathematical Physics. IBM
Journal, pages 215–234, 1967.
2. J. Etienne, M. Rastello, and E. J. Hopﬁnger. Modellingand Simulation of Powder-Snow Avalanches. C. R.
Mechanique 334, pages 545–554, 2006.
3. R. Fedkiw, J. Stam, and H. W. Jensen. Visual Simulation of Smoke. Proc. of SIGGRAPH2001, pages 390–
399, 2001.

4. Y. Gao, Li, S. C. Hu, and B. A. Barsky. Simulating
Gaseous Fluids with Low and High Speeds. Paciﬁc
Graphics 28, pages 1845–1852, 2009.
5. R.A. Gingold and J.J Monaghan. Mon. Not. R. Astron.
Soc., 181, 375.
6. J. Hong and C. Kim. Discontinuous Fluids. Proceedingsof ACM SIGGRAPH 2005, pages 915–920, 2005.
7. K. Hutter, M. Siegel, S. B. Savage, and Y. Nohguehi.
Two-Dimensional Spreading of a Granular Avalanche
down an Inclined Plane part I. Theory. Acta Mechanica
100, pages 37–68, 1992.
8. A. Kapler. Avalanche! snowy FX for XXX. SIGGRAPH
2003 Sketches and Applications, 2003, 2003.
9. S. Koshizuka and Y. Oka. Moving Particle SemiImplicit Method for Fragmentation of Incompressible
Fluid. Nuclear Science and Engineering 123, pages
421–434, 1996.
10. N. Maeno, Y. Endo, E. Akitani, S. Kobayashi, and
M. Takeuchi. Avalanche and Snowstorm (Nadare to
Fubuki). Kokonshoin, 2000 (in Japanese).
11. J. Magnaudet and I. Eames. The Motion of HighReynolds-Number Bubbles in Inhomogeneous Flows.
Annual Review of Fluid Mechanics 32, pages 659–708,
2000.
12. M. Müller, D. Charypar, and M. Gross. Particle-Based
Fluid Simulation for Interactive Applications. Proc.
of Siggraph Symposium on Computer Animation, pages
154–159, 2003.
13. A. Perla. Avalanche Release, Motion, and Impact. Dynamics of snow and ice masses, pages 397–462, 1980.
14. B. Salm and H. Gubler. Measurement and Analysisof the Motion of Dense Flow Avalanches. Annals of
Glaciology 6, pages 26–34, 1985.
15. B. Sovilla, S. Margreth, and P. Bartelt. On Snow Entrainment in Avalanche Dynamics Calculations. Cold
Regions Science and Technology 47, pages 69–79, 2007.
16. J. Stam. Real-Time Fluid Dynamics for Games. Proc.
of the Game Developer Conference, 2003.
17. M. Takeuchi. Vertical Proﬁle and Horizontal Increase ofDrift Snow Transport. Journal of Glaciology 26, pages
481–492, 1980.
18. M. Takeuchi. An Experimental Study on Tensile Fracture of Snow. Institute of Low Temperature Science,
pages 1–37, 1983.
19. G. Turk. Generating Textures on Arbitrary Surfaces Using Reaction-Diffusion. Proceedings of the 18th annual
conference on Computer graphics and interactive techniques, pages 289–298, 1991.
A Approximation on Probability of Snow Entrainment
Calculation
As shown in Section 4.4.2, we approximate that Pe is linearly proportional to ∆ t. The derivation for the approximation follows the following

9
steps. We ﬁrst consider the following two cases: 1) the simulation proceeds with large time step ∆ t, 2) the simulation proceeds with small
time step ∆ t/m in m times, where m is a positive integer. Next, assuming the snow entrainment occurs (i.e., the strain rate is greater than thethreshold), we calculate the probabilities of a single particle of the accumulated snow layer being entrained during the time ∆ t in these twocases. In order to ensure the entrainment process not to be affected by
the setting of the time step, these two probabilities should be the same.
From this observation, we ﬁnally derive the relationship between Pe
and the time step.
We ﬁrst consider the ﬁrst case. Assume that we conduct the determination of entrainment in probability P1 . Obviously, the particle isentrained during the time ∆ t in the probability P1 . Next, we consider
the second case. Assume that we conduct the determination of entrainment in a single time step ∆ t/m in probability Pm . Then, the particle isentrained during the time ∆ t in the probability Qm , given as:
Qm = Pm + (1 − Pm )Pm + (1 − Pm )2 Pm + · · · + (1 − Pm )m−1 Pm .

(17)

This equation is obtained by simply summing up the entrainment probability at each time step, and simpliﬁes to Eq. (18).
Qm = 1 − (1 − Pm )m .

(18)

When Pm is sufﬁciently small, which is the case in our simulation, Eq.
(18) can be approximated as
Qm ≈ 1 − (1 − mPm ) = mPm .

(19)

Since the probabilities of the particle being entrained during the time
∆ t in these two cases should be the same, we obtainPm ≈ P1 /m,(20)

which states that the probability Pm should be proportional to the time
step.

10

Y USUKE T SUDA received the B.S. and M.S. degrees in Computer Science from the University of Tokyo, Japan, in 2008, and 2010, respectively. His research interests center in computer graphics including visual simulation.
YONGHAO Y UE received the B.S. and M.S. degrees in Computer Science from the University of Tokyo, Japan, in 2005, and 2007, respectively. He is currently a doctoral student in the Department of Information Science and Technology at the University of Tokyo. His researchinterests center in computer graphics including photo-realistic rendering.
YOSHINORI D OBASHI is an associate professor at Hokkaido University in the graduate school of information science and technology,
Japan since 2006. His research interests center in computer graphics including realistic Image synthesis, efﬁcient rendering, sound modelingfor virtual reality applications. Dobashi received his BE, ME and Ph.D
in Engineering in 1992, 1994, and 1997, respectively, from Hiroshima
University. He worked at Hiroshima City University from 1997 to 2000
as a research associate.
T OMOYUKI N ISHITA received the B.E., M.E., and Ph.D. degrees inElectrical Engineering from Hiroshima University, Japan, in 1971,
1973, and 1985, respectively. He worked for Mazda Motor Corp. from1973 to 1979. He has been a lecturer at the Fukuyama University since
1979, then became an associate professor in 1984, and later becamea professor in 1990. He moved to the Department of Information Science of the University of Tokyo as a professor in 1998 and now is a
professor at the Department of Complexity Science and Engineering of
the University of Tokyo (since 1999). He received a Research Award
on Computer Graphics from IPSJ in 1987, Steven A. Coons Awards
from ACM SIGGRAPH in 2005, and a CG Japan Award in 2006. His
research interest is mainly in computer graphics. He is a member of
IEICE, IPSJ, ACM, and IEEE.