Visual Simulation of Compressible Snow with Friction and Cohesion
Tetsuya Takahashi1)

Issei Fujishiro2)

Tomoyuki Nishita3)

1) UEI Research / Graduate School of Science and Technology, Keio University
2) Graduate School of Science and Technology, Keio University
3) UEI Research / Faculty of Economic Information, Hiroshima Shudo University
E-mails: tetsuya.takahashi@uei.co.jp
fuji@ics.keio.ac.jp
tomoyuki.nishita@uei.co.jp

Abstract
Recent advances in physically-based simulation have made it possible to simulate various kinds of natural phenomena. However, characteristic behavior of accumulated snow while being compressed due topressure has been insuﬃciently simulated merely with a height field approximation. In this paper, we
propose a new method for simulating compressible snow by approximating snow as a set of porous snow
particles. We introduce a new parameter called durability for each porous snow particle, which represents
the rate of undamaged snow structures. The compression of accumulated snow is achieved by absorbing
the impact from solids on the snow by taking the durability of snow particles into account. In addition
to the compressibility, we incorporate friction and cohesion into our simulation framework to represent
various eﬀects of accumulated snow. For enhanced visual reality, fine scale simulations are performed asa post-processing. Several examples demonstrate the versatility of our method.

1

Introduction

Snow can be seen in various places and deeply linked toour daily lives including disasters, entertainments, andsports in cold regions. Visually appealing and attractive
snow scenes have been frequently adopted in applications such as feature films and video games. However,simulating accumulated snow still remains as a challenging issue due to its complexity. For example, theproperties of constituent particles of accumulated snow
can be easily changed by their own weight and the influence from the surrounding environment, and hence thecharacteristics of the snow drastically vary according to
the state of the particles. Because it is diﬃcult to handle such complicated conditions, three dimensional dynamics of the accumulated snow has been insuﬃciently
simulated merely with a height field approximation.
Since air volume occupies about 90% of the wholevolume of freshly-fallen snow according to Nakaya [14],newly accumulated snow collapses releasing air from inside when strongly pressed, and consequently loses theoriginal volume of the snow. Although such compression is one of the most important features of the accumulated snow, fast and simple methods for fully addressing this issue have not yet been proposed.
In this paper, therefore, we propose a new method

for simulating three dimensional and compressible accumulated snow. We introduce a new parameter called
durability for each of porous snow particles, to indicate
how much snow structures remain without being damaged, and approximate snow volume as a continuumwith these particles. In order to simulate the dynamics
of accumulated snow, we employ Fluid-Implicit-Particle(FLIP) [22] as our underlying continuum solver since
the accumulated snow partly behaves as with inviscid
fluid, viscous fluid, or granular materials. We realize thecompression of accumulated snow by reducing the eﬀectof the pressure from solids impacting on the snow taking
the durability of porous snow particles into account. We
also incorporate friction and cohesion, which are important features of accumulated snow, to improve the versatility of our method. Friction of accumulated snow issimulated by employing the Drucker-Prager yield condition and considering local pressures in the snow volume,and cohesion is achieved by the control of a durabilityrelated parameter in the yield condition. We perform
fine scale simulations as a post-processing to improvethe visual quality of accumulated snow. Figure 1 illustrates important features of our method: compression(left), friction (middle), and cohesion (right).

Figure 1: (Left) a block of snow in a gray box compressed without scattering by a solid ball moving along a
vertical path. (Middle) a ball-shaped frictional white powder dropped onto a red solid ball. (Right) a bunnyshaped white powder dropped onto the ground in a virtual box, while preserving its original shape due to thecohesion model, making several shadows.

2

Related Work

posed.

In Section 2.1, we briefly survey previous attempts torepresent snow scenes. Then in Section 2.2, we focus on 2.2 Granular Materials
simulation methods for granular materials, with which Bell et al. [5] proposed a molecular dynamics method,
our method partly shares the common features.
which is similar to Discrete Element Method (DEM), forsimulating granular materials. Although their method
is versatile and enabled two-way coupling with rigid
2.1 Snow
bodies composed of granular particles, the number ofIn the field of computer graphics, many methods for simulation particles is likely to be extremely large, andrepresenting snow have been developed with primary makes this method impractical.
focus on realistic scenes of accumulated snow. Nishita
Li and Moshell [13] proposed a height field method
et al. [16] modeled accumulated snow by making use of which approximates massive grains in order to reduce
metaballs and rendered them taking optical characteris- the computational cost. Sumner et al. [21] extended
tics which cause multiple scattering into account. An- this method for generating footprints and tracks, andother well-known method for simulating accumulated Onoue and Nishita [17] developed multi-valued heightsnow was proposed by Fearing [6], who used polygons fields and partly simulated three dimensional eﬀects.
and particles, and enabled generation of physically- Although methods using a height field approximation
plausible snow scenes by conducting recursive stabil- can eﬃciently capture general motions of granular maity tests for the particles on the polygons. Feldman terials, visually-attractive three dimensional eﬀects are
and O’Brien [7] extended the method of Fearing [6] by limited.
considering air stream and produced wind-driven accuTo suﬃciently and eﬃciently simulate the three dimulated snow scenes. A level set approach was pro- mensional eﬀects of granular materials, a continuum
posed by Hinks and Museth [10] for smoothly accumu- approach was introduced by Zhu and Bridson [22], who
lated snow. Festenberg and Gumhold [8] employed a simulated sand using FLIP. Narain et al [15] simulated
statistical snow deposition model using diﬀusion equa- cohesionless granular materials improving the method
tions, and successfully generated snow scenes. Tsuda in [22]. Alduán and Otaduy [1] and Ihmsen et al.
et al. [20] simulated avalanche with snow smoke. Taka- [11] generated motions of granular materials in thehashi and Fujishiro [19] enabled the simulation of adhe- Smoothed Particle Hydrodynamics (SPH) framework.
sive snow interacting with solid objects. Stomakhin etal. [18] proposed a new method which uses Materials
Point Method (MPM) and simulated the dynamics of 3
Simulation Framework
accumulated snow in a sophisticated manner.
Though there are many well-known methods for sim- We assume that accumulated snow is a continuum sinceulating accumulated snow, fast and simple methods for individual snow particles are quite tiny, and thus theirsuﬃciently addressing the compressibility of three di- very detailed dynamics is visually-ignorable. We hereinmensional accumulated snow have not yet been pro- employ FLIP [22] as our continuum solver for generat-Figure 2: An overview of our method, where gray lines represent a grid, a black box is a solid body, blue particles
are porous snow ones, and red particles are fine scale ones. These objects are translucently drawn when they
are irrelevant to the simulation in each step. Our simulation has a loop consisting of two phases, each of which
has two steps. In the particle phase, particle positions are updated while the fine scale simulation is performedat the first step, and velocity and durability of particles are transferred onto the grid at the second step. In thegrid phase, computations of snow motions are conducted by updating pressure, friction, and cohesion at the firststep, and grid velocity and durability are interpolated to the particles with FLIP at the second step.

ing snow motions because the solver is robust regardless We employ an SPH-like interpolation as in [4, 2], andof particle distributions unlike SPH and DEM, and this A and ρ at position x are transferred as follows:feature is necessary for our goal. Figure 2 illustrates an∑i mi Ai Wsharp (xi − x, h)overview of our method. Our method consists of two,(2)
A(x) = ∑phases (particle phase and grid phase), each of which
i mi Wsharp (xi − x, h)∑
has two steps, and the phases are executed alternately
ρ(x) =
mi Wsmooth (xi − x, h),and repeatedly.
i

3.1

Fundamentals

The motion of a continuum is represented with density
ρ and velocity u. If we denote the Lagrangian derivativeby D/Dt = ∂/∂t + (u · ∇), Navier-Stokes equations are
given asρDu
= −∇p + ∇ · s + fext ,Dt(1)

where i denotes the index of a particle within the effective radius h. With a vector r, we use the following
kernels proposed in [2] because of their simplicity andeﬃcacy:{ 2
h /||r||2 − 1
0 ≤ ||r|| ≤ h,
Wsharp (r, h) =0otherwise,{2
2
1 − ||r|| /h
0 ≤ ||r|| ≤ h,Wsmooth (r, h) =0otherwise.

where p denotes pressure, s deviatoric or frictional We employ the ghost fluid method [9] and variational
stress, and fext external forces. We obtain the dynamics framework [3] for accurate simulations.
of the continuum by solving Eq. (1) with Dirichlet (forfree surfaces) and Neumann (for solid objects) bound- 3.3 Fine Scale Simulation
ary conditions, and the incompressibility constraint
In our method, since simulation particles represent∇ · u = 0.
clamp of matter or distributions of a continuum insteadof individual snow particles, we do not directly renderthese simulation particles. For particles used in render3.2 FLIP Solver
ing and enhanced visual reality, we perform a fine scale
In FLIP [22], particle properties are transferred onto simulation as a post-processing [11].
the grid. Let A denote an arbitrary physical quantity
In the fine scale simulation, we ignore interactions
except density, m mass of a particle, and W a kernel. among particles for eﬃciency. Fine scale particles aresŽŝĚƐ
ƵƌĂďŝůŝƚǇ݀
^ŶŽǁ
^ŶŽǁ

DŝĐƌŽƐĐŽƉŝĐƐŶŽǁDĂĐƌŽƐĐŽƉŝĐƐŶŽǁ

Figure 3: Definition of a porous snow particle.
Figure 4: Snow balls with diﬀerent durabilities interacting with moving cuboids.
governed by the dynamics of the continuum at highdensity points, and can freely disperse due to external
forces in low density areas. Hence, we compute the pressure from the solid object due to the voids acting as
velocity of fine scale particles v by taking external forces buﬀers. Consequently, the velocities of the snow volume
and the solid object are not the same, and the non-zerofext into account:relative velocities make the snow get compressed.
v = αu + (1 − α)(u + fext ∆t/m),
To realize the behavior of compressible accumulated{snow, we transfer the durability of particles to the gridρ/ρ∗ρ ≤ ρ∗ ,α=using Eq. (2), and for each grid cell, compute the ve1
otherwise,
locity transfer factor b (0 ≤ b ≤ 1), which controls thewhere ∆t denotes a time step and ρ∗ a threshold density. extent of influence of solid velocity to the continuum
motions by

4

Compressibility

Snow has compressibility, due to air voids in its volume, which is not exhibited by inviscid fluid, viscous
fluid, nor granular materials. Section 4.1 describes how
to handle these voids in our model, and Section 4.2
explains the eﬀects of colliding objects to realize thecompression of accumulated snow.

4.1

Air Voids

Most part of a snow volume is occupied by air since
it has myriad air voids in a sampled volume unlike
other continua. Accumulated snow, therefore, possesses
a characteristic trait that it loses the volume releasing
air stored inside when being pressed. To simulate thisphenomenon, we utilize porous snow particles, whichwere inspired by the method in [12]. As illustrated in
Figure 3, we handle snow particles in a macroscopic
scale by parameterizing the detailed air voids and microscopic snow structures as a scalar variable, durabilityd (0 ≤ d ≤ 1), which indicates how much the original snow structures remain without damage. Note thatdurability d = 1 means that snow structures of a particle are intact.

4.2

Eﬀects of Colliding Objects

When a solid object collides with accumulated snow,
the snow absorbs the impact owing to its voids whilesnow structures are damaged. Therefore, the velocity ofthe snow volume is not likely to be fully aﬀected by theb = b0 − kb (d0 − d),where b0 denotes an initial velocity transfer factor,
kb (0 ≤ kb ) a coeﬃcient of velocity transfer factor, andd0 an initial durability. Then, a Poisson equation ismodified as∆t 2∇ p = ∇ · (u − bus ),ρ
where us denotes the solid velocity. Figure 4 illustrates
the eﬀect of durability with an example of two snowballs interacting with two moving cuboids, respectively.
The left snow ball with high durability absorbs the impact due to the left cuboid, whereas the right one withlow durability is strongly scattered.
Since when accumulated snow pressed by an object,
the structures of the snow are damaged. We assume
that the entire energy generated by the collision is used
for the destruction of snow structures ignoring the dissipation of the energy as sound and heat for simplicity,and then reduce the grid durability d byd ← d − kq p,(3)

where kq denotes a durability change coeﬃcient. Sincedamaged structures never recover, the durability reduces monotonously. After all the durability values on
the grid are updated with Eq. (3), they are interpolated with FLIP [22]. Note that Particle-In-Cell (PIC)[22] is inapplicable to the interpolation due to the fast
numerical dissipation.

5

Friction and Cohesion

As in the simulations for granular materials, a yield
condition is responsible for determining the dynamics
of accumulated snow. We employ the Drucker-Prager
yield condition with a cohesion term to check whether
s is large enough to make the snow plastically flow:||s||F ≤√3αp + c,(a) µ0 = 0.0, kf = 0.00(b) µ0 = 0.0, kf = 0.02√
where ||s||F = Σs2ij denotes the Frobenius norm of s,√
α = 2/3 sinθ, where θ is an angle of repose, and c acohesion criterion. If the yield condition does not holdin a cell, the cell is labeled as plastic, otherwise rigid.
If the label of cells is plastic, the continuum in thecells plastically flows dissipating its energy owing to
the friction. To simulate such behavior, we formulate
implicit systems for friction, which are similar to theimplicit systems for viscosity proposed by Batty and(c) µ0 = 0.0, kf = 0.06(d) µ0 = 0.0, kf = 0.15
Bridson [4]. Since friction forces become stronger aspressure gets higher, we compute a friction coeﬃcientµ byFigure 5: A comparison of friction eﬀects with diﬀerent
parameter values.
µ = µ0 + kf p,
where µ0 denotes an initial friction coeﬃcient, and kf
a friction change coeﬃcient. Figure 5 demonstrates the
Table 1: Statistics of all the scenes.
eﬀect of our friction model using ball-shaped white powders with diﬀerent parameter values, dropped onto a red
Figure Grid
ParticlesTime/solid ball.
number size Simulation/fine scale simulation step (s)In order to handle cells labeled as rigid, we follow
4503
28.3k/1.7M
2.1
5
503
14.2k/0.5M
0.6
the method in [22] for cohesive motions. First, we con3
66064.8k/1.9M
3.3
nect rigid-labeled cells to construct rigid clusters, and7403
128.3k/0
0.8
compute the linear and angular velocity of the clusters
8503
224.1k/4.0M
1.4
preserving their linear and angular momentum. Then,
11503
33.0k/2.0M
1.1
we update the velocity of rigid-labeled cells with the linear and angular velocity of the rigid clusters. Accumulated snow generally shows diﬀerent cohesion depending
on positions in the snow volume. For example, if thestructures of accumulated snow are damaged, the snow 6
Results and Discussions
shows higher cohesion, while freshly-fallen and undamaged snow does not exhibit an adhesive property. Bytaking these facts into account, we compute the cohe- We implemented our method in C++ and parallelized itsion criterion c with the durability d byemploying OpenMP 2.0. All the scenes were executedon a PC with Intel Core i7 3.50 GHz and RAM 16.0
c = c0 + kc (d0 − d),GB. We used six threads in total and obtained the accelerated simulation performance by a factor from two
where c0 denotes an initial cohesion criterion, and kc a to four. All the results were rendered with a raytracer
cohesion change coeﬃcient. Figure 6 illustrates the dif- POV-Ray 3.7. We rendered particles used in the fineference of a dropped bunny with diﬀerent conditions. scale simulation and needed about 90 seconds per frame
In this scene, while the bunny with little cohesion scat- on average for oﬀ-line rendering of Figure 1 (left). We
ters as with fluid, the highly cohesive bunnies preserve tabulate the simulation conditions and performances of
all the results in Table 1.
their original shapes.

(a) c0 = 0.0, kc = 0.0(b) c0 = 20.0, kc = 0.0(a) Snow(b) Granular materials(c) c0 = 40.0, kc = 0.0(d) c0 = 60.0, kc = 0.0(c) Viscous fluid(d) Fluid

Figure 6: A comparison of cohesion eﬀects with diﬀer- Figure 7: Motion comparisons of continua interactingwith a solid cube moving along a vertical path, whereent parameter values.
simulation particles are rendered for visualization.

6.1

Results

Figure 7 compares snow (a), granular materials (b),viscous fluid (c), and fluid (d) with particle representations, where particles are compressed by the cube
moving along a vertical path. While the snow block
is compressed due to the pressure from the cube without scattering and flowing in Figure 7 (a), the granularmaterials are scattered by the cube in Figure 7 (b). In
Figure 7 (c), the viscous material is pushed by the cubeand rises avoiding the cube, and the fluid flows due tothe pressure from the cube in Figure 7 (d).
Figure 8 illustrates a snow block being compressedby the solid ball, and Figure 9 provides side views of
the scene (Figure 8) in terms of density to elucidate
the snow compression, where the front side is clippedfor visualization, and particles are colored according totheir density (low density: white and high one: red).
We can observe that the block of snow is compressed
without flowing and scattering due to our model of thecompressibility, friction, and cohesion, and that many
particles are gathered under the solid due to the compression. Figure 10 plots a sequence of timings to generate the results in Figure 8. On average, the computational time consumed for the simulation is divided asfollows: 7% for compressibility, 14% for pressure, 25%

for friction, 1% for cohesion, 35% for fine scale simulation, and 18% for the rest. Most of the simulation time
is used for the computations of friction and fine scale
simulation, and additional cost for the compressibility
in our method is relatively low.
Figure 11 (left) demonstrates a snow ball dropped on
the solid cube, where the snow stays on the box thanks
to the friction and cohesion forces. A corresponding
simulation particle representation is shown in Figure
11 (right).

6.2

Discussions

Unlike FLIP [22] that we used, it would be diﬃcult
to capture the behavior of compressible snow using aparticle-based method such as SPH or DEM. It is because particle distributions in the particle-based methods can have much impact on the motions of the particles, and the particles could cause extremely strong
repulsion forces to resolve penetrations and collisionsof the particles when compressed by excessively heavy
objects. Although adopting very small time steps or
merge/split operations of particles would alleviate theproblem, a significant increase of computational cost
would be anticipated. On the other hand, FLIP [22] is900
Compressibility
Pressure
Friction
Cohesion
Fine scale simulation

800

Simulation time (ms)

700
600
500
400
300
200
100
0
0

50

100
Simulation steps

150

200

Figure 10: A timing profile of the snow compression in
Figure 8.

Figure 8: A block of snow compressed by a solid ballmoving along a vertical path. Simulation steps are 0,
18, 36, and 132.

Figure 11: (Left) a snow ball dropped onto a box.
(Right) the corresponding simulation particle representation.

Figure 9: Side views of the snow compression showed in
Figure 8, where the front side is clipped for visualization, and particles are colored according to their density(low density: white and high one: red).

fairly tolerable to biased distributions of particles since
FLIP [22] does not need to resolve particle penetrations
and collisions.

7

Conclusions and Future Work

We have presented a new method for simulating three
dimensional and compressible accumulated snow. Our
method utilizes porous snow particles, which have durability, to allow for the compression of accumulated snow
by the control of the impact of solids on the snow. Additionally, we modeled frictional motions of accumulated
snow formulating implicit systems, and realized cohesion of snow by taking the durability into account.
In our method, however, accumulated snow cannot

freely disperse and gradually flows as with viscous fluid
since we assume that a continuum which approximates
the snow is incompressible. Narain et al. [15] modified the incompressible FLIP and proposed the unilateral incompressibility constraint to simulate freelyflowing granular materials by solving the linear complementarity problem. Stomakhin et al. [18] utilized
MPM to simulate the dynamics of accumulated snow,
and produced scattering snow particles. Although our
method focused on simple models for fast simulations,these complicated methods would be promising to realize dispersing particles in our framework. In addition,
we would like to simulate the phase transition of accumulated snow caused by pressure melting, and rendersnow particles taking optical properties into account.

References[1] Alduán, I. and Otaduy, M. A.: “SPH Granular
Flow with Friction and Cohesion,” in Proceedingsof the 2011 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, pp. 25–32, 2011.
[2] Ando, R., Thürey, N., and Turuno, R.: “Pre-serving Fluid Sheets with Adaptively Sampled [13] Li, X. and Moshell, J. M.: “Modeling Soil:
Realtime Dynamics Models for Soil Slippage andAnisotropic Particles,” IEEE Transactions on ViManipulation,” in Proceedings of SIGGRAPH, pp.
sualization and Computer Graphics, Vol. 18, No.
8, pp. 1202–1214, 2012.
361–368, 1993.
[3] Batty, C., Bertails, F., and Bridson, R.: “A [14] Nakaya, U.: Snow Crystals: Natural and Artificial,
Harvard University Press, 1954.
Fast Variational Framework for Accurate SolidFluid Coupling,” ACM Transactions on Graphics, [15] Narain, R., Golas, A., and Lin, M. C.: “FreeVol. 26, No. 3, Article 100, 2007.
Flowing Granular Materials with Two-Way Solid
Coupling,” ACM Transactions on Graphics, Vol.
[4] Batty, C. and Bridson, R.: “Accurate Viscous
29, No. 6, Article 173, 2010.
Free Surfaces for Buckling, Coiling, and Rotating
Liquids,” in Proceedings of the 2008 ACM SIG- [16] Nishita, T., Iwasaki, H., Dobashi, Y., and NakaGRAPH/Eurographics Symposium on Computermae, E.: “A Modeling and Rendering Method forAnimation, pp. 219–228, 2008.
Snow by Using Metaballs,” Computer Graphics Forum, Vol. 16, No. 3, pp. 357–364, 1997.
[5] Bell,N., Yu,Y., and Mucha, P.
J.:“Particle-Based Simulation of Granular Mate- [17] Onoue, K. and Nishita, T.: “Virtual Sandbox,” inProceedings of IEEE 2003 Pacific Conference on
rials,” in Proceedings of the 2005 ACM SIGComputer Graphics and Applications, pp. 252–259,
GRAPH/Eurographics Symposium on Computer
2003.
Animation, pp. 77–86, 2005.
[6] Fearing, P.: “Computer Modelling of Fallen [18] Stomakhin, A., Schroeder, C., Chai, L., Teran, J.,
and Selle, Andrew.: “A Material Point Method for
Snow,” in SIGGRAPH ’00: Proceedings of the 27th
Snow Simulation,” ACM Transactions on Graphics
annual conference on Computer Graphics and InVol. 32, No. 4, pp. 102:1–102:10, 2013.
teractive Techniques, pp. 37–46, 2000.
[7] Feldman, B. E and O’Brien, J. F.: “Modelingthe Accumulation of Wind-Driven Snow,” in SIGGRAPH ’02: ACM SIGGRAPH 2002 Conference
Abstracts and Applications, pp. 218–218, 2002.
[8] Festenberg, N. V. and Gumhold, S.: “DiﬀusionBased Snow Cover Generation,” Computer Graphics Forum, Vol. 30, No. 6, pp. 1837–1849, 2011.

[19] Takahashi, T. and Fujishiro, I.: “Particle-Based
Simulation of Snow Trampling Taking Sintering Effect into Account,” in SIGGRAPH 2012 Posters,
Article No. 7, 2012.
[20] Tsuda, Y., Yue, Y., Dobashi, Y., and Nishita, T.:“Visual Simulation of Mixed-Motion Avalancheswith Interactions Between Snow Layers,” The Visual Computer, Vol. 26, No. 6–8, pp. 883–891,
2010.

[9] Gibou, F., Fedkiw, R., Cheng, L.-T., and Kang,
M.: “A Second-Order-Accurate Symmetric Dis- [21] Sumner, R. W., O’Brien, J. F., and Hodgins, J.
K.: “Animating Sand, Mud, and Snow,” Computer
cretization of the Poisson Equation on Irregular
Graphics Forum, Vol. 18, No. 1, pp. 17–26, 1999.
Domains,” Journal of Computational Physics, Vol.
176, No. 1, pp. 205–227, 2002.
[22] Zhu, Y. and Bridson, R.: “Animating Sand as a
Fluid,” ACM Transactions on Graphics, Vol. 24,[10] Hinks, T. and Museth, K.: “Wind-Driven Snow
No. 3, pp. 965–972, 2005.
Buildup Using a Level Set Approach,” in Proceedings of Eurographics Ireland Workshop Series, Vol.
9, pp. 19–26, 2009.
[11] Ihmsen, M., Wahl, A., and Teschner, M.: “HighResolution Simulation of Granular Material with
SPH,” in Proceedings of the VRIPHYS, pp. 53–60,
2012.
[12] Lenaerts, T., Adams, B., and Dutré, P.: “Porous
Flow in Particle-Based Fluid Simulations,” ACM
Transactions on Graphics, Vol. 27, No. 3, Article
49, 2008.