Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2010)
M. Otaduy and Z. Popovic (Editors)

Real-time Simulation of Large Bodies of Water
with Small Scale Details
Nuttapong Chentanez and Matthias MüllerNVIDIA Corporation

Abstract
We present a hybrid water simulation method that combines grid based and particles based approaches. Our specialized shallow water solver can handle arbitrary underlying terrain slopes, arbitrary water depth and supports
wet-dry regions tracking. To treat open water scenes we introduce a method for handling non-reflecting boundary
conditions. Regions of liquid that cannot be represented by the height field including breaking waves, water falls
and splashing due to rigid and soft bodies interaction are automatically turned into spray, splash and foam particles. The particles are treated as simple non-interacting point masses and they exchange mass and momentum
with the height field fluid. We also present a method for procedurally adding small scale waves that are advectedwith the water flow. We demonstrate the effectiveness of our method in various test scene including a large flowing
river along a valley with beaches, big rocks, steep cliffs and waterfalls. Both the grid and the particles simulations
are implemented in CUDA. We achieve real-time performance on modern GPUs in all the examples.
Categories and Subject Descriptors (according to ACM CCS): Computer Graphics [I.3.5]: Computational Geometry
and Object Modeling, Physically Based Modeling—Computer Graphics [I.3.7]: Three-Dimensional Graphics andRealism, Animation—Simulation and Modeling [I.6.8]: Type of Simulation, Animation—
Keywords: natural phenomena, physically based animation

1. Introduction
In recent years, physically based effects have become more
and more popular in computer games. Features that used tobe hand-animated are now driven by simulators. The mostimportant examples are rigid bodies and particles. However,in most cases, large bodies of water such as rivers, lakes oroceans are still animated using procedural approaches. The
reason for this is the large simulation domain in these scenarios. To get details on the surface of a lake using a full 3D
fluid simulation would require millions of grid cells in the
Eulerian case or millions of particles if an SPH approach isused. For real-time applications such as computer games, aviable simulation approach is to reduce the problem from 3D
to 2D by treating the body of water as a 2D height field. This
reduces simulation time substantially but comes at a price.
Phenomena like waterfalls or overturning waves cannot be
represented by a 2D height field. We solve this problem byturning height field water into particles automatically basedon certain criteria that depend on the shape of the surface.
c The Eurographics Association 2010.

These particles are turned back into height field water when
they hit the water surface.
A simple and quite popular technique to simulate water surfaces is to solve the 2D wave equation. This approach works
well for local phenomena like puddles or waves around aboat. The main drawback of this mathematical model, however, is that it only works with a vertical velocity field.
Therefore, effects that are based on horizontal motion suchas whirlpools or the flow of a river cannot be treated correctly. This is the reason we decided to use the shallow waterequations for our simulations. The presence of a horizontal
velocity field allows us to correctly advect floating objects
and foam particles on the surface.
In the following sections we will show how various techniques can be put together to create a comprehensive system for simulating large bodies of water with small scaledetail in real-time. This paper is not a system paper in the
traditional sense though because we did not simply combine

Nuttapong Chentanez & Matthias Müller / Real-time Simulation of Large Bodies of Waterwith Small Scale Details

well known components but came up with various new approaches to solve the problems that arise when attempting tosimulate a realistic water scenario with all its large and small
scale features.

Staadt [CS09]. The main drawback of methods based on the
wave equation or the pipe model is that vortices which areresponsible for many interesting water phenomena such as
swirling foam and whirlpools are not present in the model.

The main contributions of this paper are:

The Shallow Water Equations (SWE) can capture these phenomena because in addition to the height field, they describe the evolution of a 2D velocity field normal to the water columns. The SWE were introduced to computer graphics in [LvdP02]. The approach is based on the idea of assuming linear vertical pressure profiles in the 3D NavierStokes equations [CL95]. Hagen et al. [HHL∗ 05] used afinite volume method for solving the SWE on the GPU tosimulate water flowing on irregular terrains. As height field
fluid approaches neither the wave equation nor the SWE
can handle breaking waves. [TMFSG07] proposed a method
to simulate breaking waves by automatically generating triangle mesh patches and evolving them. Their method cannot be parallelized easily though. [WMT07] added surface
tension forces to the SWE to simulate water flow on arbitrary surfaces. Later, [ATBG08] used an implicit Newmark integration scheme to reduce numerical dissipation in
the velocity advection step. To simplify the solver, [LO07]
ignored the divergence term and used a collocated grid.
Hess et al. [Hes07] integrated the SWE explicitly and described methods for two way water - rigid body coupling
and wet-dry region tracking. The grid based portion of our
proposed method also employs the SWE with several novel
enhancements. Other approaches to model larger bodies ofwater such as rivers include the use of stream functions in[YNBH09] and the use of the 2D Navier-Strokes equations
coupled loosely with the pipe model in [BAB09]. They advect particles that store texture coordinates for adding small
scale detail to the water.

• Various improvements to a basic shallow water solver thatenhance stability and allow the handling of arbitrary terrain slopes and arbitrary water depth.
• Criteria for automatically detecting surface discontinuities such as waterfalls and a method to turn the liquidvolume in such locations into particles that carry mass andmomentum of the height field across the discontinuity.
• A technique to automatically create and advect small scale
FFT-based ripples on the water surface.
• A new and robust criterion to detect breaking waves anda method to turn height field fluid into particles creating
spray, splash and foam.
• A way to absorb waves at the boundary of the simulation
domain for handling open water scenes.
1.1. Related Work
Procedural methods for generating water animations havebeen used since the early days of computer graphics [FR86,
Pea86]. There are several papers that describe how to use thespectral method for modeling ocean waves based on the Fast
Fourier Transformation (FFT) in order to create animationsof large bodies of water [MWM87, Tes99, TDG00, HNC02].
Although these methods are well suited to generate high resolution and large scale water animations, they cannot model
the interaction with solid objects easily and are unable togenerate vortices. With the exception of [Pea86] the separation of liquid from the main body of water has not been
considered in these papers.
In computer graphics, [KM90] and [Tes99] were among thefirst to simulate water surfaces by solving the wave equation with internal boundaries on a 2D height field. Many
researchers have used a pipe model where adjacent cells
are connected by pipes that allow water to flow across cellboundaries. O’Brien and Hodgins [OH95] extended the pipe
model by adding splash particles and modeling rigid bodyinteraction with the water. Mould and Yang [MY97] later
added bubbles and modeled droplets. Holmberg and Wünsche [HW04] used a weir model to allow the simulation ofwaterfalls. The method was later accelerated by a GPU implementation in [MFC06]. The pipe model was also used tosimulate hydraulic erosion of terrains [vBBK08]. The authors produced impressive results by accelerating the simulation on GPUs. Recently, [YHK07] introduced a novel approach to approximate the solution of the wave equation bystoring wave trains on 2D particles. Multiple shifting grids
are used for simulating open water scenes by Cords andThere are many works on the full 3D simulation of water
which are beyond the scope of this paper. We wholeheartedly refer the reader to the book by Bridson [Bri08] forgrid based simulation, and to a paper by Solenthaler and Pajarola [SP09] and references therein for particles based simulation.
2. Methods
We will state all the steps explicitly in full algorithmic detail.
The goal is to provide enough information for the reproduction of our results. To get an overview of our method, thereader can skip the formulas and concentrate on the descriptions in the text.
Algorithm 1 Main loop (one time step)
1. Height field fluid simulation (Section 2.1)
2. Solids simulation (any standard simulator)
3. Two-way coupling of height field and solids(Section 2.3)
4. Particles generation and simulation (Section 2.4)
5. Rendering (Section 2.5)c The Eurographics Association 2010.

Nuttapong Chentanez & Matthias Müller / Real-time Simulation of Large Bodies of Waterwith Small Scale Details

Our main loop is summarized in Algorithm 1. We first simulate the height field fluid. Gravity, external forces and boundary conditions are taken into account as discussed in Section 2.1. Then, we simulate the solids such as rigid and soft
bodies using a standard simulation package. We then couple
the solids and the height field fluid as discussed in Section
2.3. Next, forces that the fluid exerts on the solids such as
buoyancy, drag and lift are computed and applied. We thenmodify the fluid’s height field and velocity field to take themovement of the solids into account. After that, we generateparticles to replace regions of the liquid that the height fieldcannot resolve well. We then simulate the movement of theparticles and eliminate them if they fall back into the height
field (see Section 2.4). Finally, we render the height fieldwith an additional displacement map that represents small
waves with wave lengths below the resolution of the simulation grid. The splash particles are rendered with a screenspace technique as explained in Section 2.5

2.1. Height Field Fluid Simulation
We employ the Shallow Water Equations (SWE) which simplify the full Navier-Stokes equations to a 2D height field
representation of the liquid surface. We assume gravity to
act along the y-axis. Hence the plane for the 2D simulationis the x-z plane. In the following equations, h is the depth of
the water, H is the y-coordinate of the terrain on the bottom,
η = H + h is the y-coordinate of the water surface, v = (u, w)
is the horizontal velocity of the fluid, g is gravity and aext is
an external acceleration. The shallow water equations can be
written asDh= −h∇ · v,
Dt
Dv
= −g∇η + aext ,
Dt

(1)
(2)

where D is the material derivative operator. The equations
describe conservation of mass and conservation of momentum. A detailed derivation of the SWE can be found in[Bri05]. Several authors [LvdP02], [WMT07], [ATBG08]
solved the SWE using implicit methods. This guarantees unconditional stability. However, they either require an iterative method which incurs a significant run time cost or needto pre-factorize the matrix which means boundary conditions
cannot be changed. We therefore use explicit integration andhandle the terms that can lead to instability specially.
We discretize the simulation domain with a staggered gridwhere the heights hi, j and Hi, j are stored at the cell centers and the velocities components ui+ 1 , j , wi, j+ 1 on faces
22
following [Bri05]. We employ a time-splitting technique by
first solving the self advection of the velocity field and thenintegrating the height field and velocity field forward in time.
Throughout this paper we use meter (m) for distance and second (s) for time as our units. The grid spacing and time steps
are denoted by ∆x and ∆t, respectively.
c The Eurographics Association 2010.

2.1.1. Velocity Advection
We solve the advection of ui+ 1 , j , vi, j+ 1 using an un2
2
conditionally stable modified MacCormack method as in∗[SFK 08] and fall back to the semi-Lagrangian method if theresulting value is not within the bounds of the velocity values
used for bilinear interpolation of the first semi-Lagrangian
sub-step as suggested in [SFK∗ 08].
2.1.2. Height Integration
Equation 1 can be re-written as∂h= −∇ · (hv),∂t(3)

and discretized it to get(h̄u)i+ 1 , j−(h̄u)i− 1 , j (h̄w)i, j+ 1 −(h̄w)i, j− 1∂hi, j
2
2
2
2
=−(
+
),(4)
∂t
∆x
∆x
where h̄ is h evaluated in the upwind direction(
hi+1, j if ui+ 1 , j ≤ 02h̄i+ 1 , j =2hi, jif ui+ 1 , j > 0,

(5)

2

and
h̄i, j+ 1 =2(
hi, j+1

if wi, j+ 1 ≤ 0hi, jif wi, j+ 1 > 0.

2

(6)

2

We integrate the height explicitly using
hi, j +=

∂hi, j∆t.
∂t(7)This guarantees mass conservation for the height integration
step. We also found that it yields a more stable simulation
h +h
h +h
than the choice of h̄i+ 1 , j = i, j 2 i+1, j and h̄i, j+ 1 = i, j 2 i, j+1
2
2
of [Hes07] in practice. The height integration step will be
modified to take the discontinuity due to waterfalls into account in Section 2.4.3.
2.1.3. Velocity Integration
The face velocities are updated by taking the gradient of thewater height into account as follows:−g(ηi+1, j − ηi, j ) + aextx )∆t,∆x−g
wi, j+ 1 += ((ηi, j+1 − ηi, j ) + aextz )∆t.
2∆xThis step will also be modified in Section 2.4.3.
ui+ 1 , j += (

(8)

2

(9)

2.1.4. Boundary Conditions
For a face that is marked by the user as reflective such as
walls or permanent static obstacles, the velocity value is always set to 0 at the end of every time step and is not updated
in the advection step. Moreover, we treat a face (i + 12 , j) as
reflective if either of the following is true:

Nuttapong Chentanez & Matthias Müller / Real-time Simulation of Large Bodies of Waterwith Small Scale Details• hi, j ≤ ε and Hi, j > ηi+1, j or
• hi+1, j ≤ ε and Hi+1, j > ηi, j ,where ε > 0 is a small constant. We use 10−4 ∆x in all examples. A similar condition applies to a face (i, j + 12 ). Conceptually, the water level in a wet cell needs to be higher than
the terrain level in the neighboring dry cell before the flow
starts. Otherwise, the face behaves as a wall.
To simulate open water scenes, borders that absorb all incoming waves are needed. We adapt the Perfectly Matched
Layers (PMLs) method [SM09], [Joh08] with a few modifications to take the rest water level into account and to improve stability. The idea is to define a damping region near
the border of the domain and make waves within that region decay fast enough such that their amplitudes are negligible when they hit the border. To damp out the waves thatmove along the x-axis we evolve a field φ, defined only in thedamping region. The field is initialized to 0 and then used to
update h and u ashi, j += (−σi, j (hi, j − hrest ) + φi, j )∆t,1
ui+ 1 , j += − (σi+1, j + σi, j )ui+ 1 , j ∆t,2
22(wi, j+ 1 − wi, j− 1 )22∆t,
φi, j += −λupdate σi, j∆x
φi, j ×= λdecay ,Figure 1: A domain with a reflecting boundary condition onleft and right, and absorbing region on the top and bottom.
The absorbing regions are highlighted in pink.

(10)
(11)
(12)
(13)

where the constant σi, j controls how fast the wave’s amplitudes decay, hrest is the depth of the water at rest and0 < λdecay < 1 and 0 < λupdate < 1 are constants to help
improve the stability of the explicit Euler time integration.
As suggested in [Joh08], σi, j should increase quadraticallyor cubically from the inner edge of the damping region to
the domain border. Large σi, j will result in a faster decay
and allow narrower damping regions but may introduce numerical reflection when a wave is about to enter the region. We use λdecay = 0.9 and λupdate = 0.1 for all the examples. Waves moving in another directions are damped by theabove method with an attenuation factor proportional to thecosine of the angle of the direction with the x-axis [Joh08].
To sufficiently damp the waves moving in all directions, asimilar formulation is used for the z-axis, see Appendix 3.1.
In all of our examples that require it, the width of the damping region is 10 cells. An example scene showing the PML
is shown in Figure 1. The top and the bottom borders areabsorbing z-moving waves while the left and the right borders are reflective. Notice how waves reflect only from theleft and right borders. The technique is used in the boat riding scene in Figure 8 and in the open ocean scene shown in
Figure 2.
2.1.5. Stability Enhancements
We propose several measures to improve the stability of thesimulation without losing much visual quality. Due to numerical error, hi, j can become smaller than zero which is not

Figure 2: An open ocean scene using an absorbing layer
around the SWE domain which is surrounded by an FFTbased ocean simulation.
physical and can cause stability issues. Therefore, we clamphi, j to always be ≥ 0.
We also clamp the magnitudes of ui+ 1 , j and vi, j+ 1 to be less22than α ∆x∆t , where we use α = 0.5 in all examples. This places
a limit on the wave speed within the simulation, but greatly
improves the stability for scenes with violent water.
Additionally, we limit the water depth used for the height
integration by replacing h̄ in Equation 4 by h̄ − hadj ,
h̄

+h̄

+h̄

+h̄

where hadj = max(0, i+1, j i−1, j 4 i, j+1 i, j−1 − havgmax ) and∆x
havgmax = β g∆t. We use β = 2.0 in all examples. This places
a limit on the wave lengths and amplitudes in deep water regions, but improves the stability of the simulation while still
produce plausible results.
2.2. Overshooting Reduction
When a wave from deeper waters enters a shallower region,
its amplitude increases and the wavelength decreases. Eventually, the wave lengths become smaller than ∆x and cannot
be resolved by the grid anymore but become numerical ripples instead. The use of explicit integration makes the situation worse as it tends to amplify the edge of the first wave
c The Eurographics Association 2010.

Nuttapong Chentanez & Matthias Müller / Real-time Simulation of Large Bodies of Waterwith Small Scale Details

v are obtained by barycentric interpolation from the corresponding original triangle. Let n = [nx , ny , nz ]T be the normal vector of the triangle. Throughout the section, ŷ denotes
the unit vector [0, 1, 0]T .
2.3.1. Fluid to Solids

Figure 3: a) Overshooting artifacts on top of breaking
waves. b) Artifacts reduced by our technique.

+-Figure 4: a) Division of triangles into similar triangles b)
Estimating the submerged volume by adding/subtracting theprism between downward/upward facing triangles and thewater surface.

We model the forces that the fluid exerts on a solid by taking three major components into account namely buoyancy,drag and lift forces. Buoyancy is an upward pointing force
proportional to the weight of the displaced fluid. It can becomputed as −gρV , where V is the volume of the object below the water surface. V is estimated simply by adding andsubtracting the volumes of the prisms formed by projecting
the triangle along the y-axis to the water surface. Volume isadded for the downward facing triangles and subtracted forthe upward facing triangles, as shown in Figure 4b.
The contribution of a sub triangle to the buoyancy force is
hence
(
0
if η(p) < pybouyancy(14)f=gρA(η(p) − py )ny ŷ otherwise,
where η(p) is the water level at p, evaluated by a bi-linear
interpolation.
The drag and lift forces are computed as in [YHK07],

front that enters the shallower water region. This manifests
itself as overshooting artifacts in the simulation result, asshown in Figure 3a. To reduce the problem, we detect edgesof these waves and reduce their magnitude by checking• if ηi, j − ηi−1, j > λedge and ηi, j > ηi+1, j , thenhi, j += αedge (max(0, 21 (hi, j + hi+1, j )) − hi, j ).
• if ηi, j − ηi+1, j > λedge and ηi, j > ηi−1, j , thenhi, j += αedge (max(0, 12 (hi, j + hi−1, j )) − hi, j ).
We use λedge = 2∆x and 0.1 < αedge < 0.5 in our examples.
A similar fix is used to waves moving along the z-axis. This
reduces the overshooting for waves traveling in other directions as well, albeit less effective. Notice how the overshooting artifact is reduced in Figure 3b.
2.3. Two-Way Coupling of the Height Field with Solids
In this section we describe how we handle the coupling of
rigid and soft bodies with the height field fluid. For all the triangles of the rigid bodies and the soft bodies we recursively
divide each triangle into smaller similar sub-triangles until
the area of these sub-triangles falls below κ∆x2 (see Figure4a). With decreasing κ the force computation becomes moreaccurate but also more expensive. We use κ = 1 in all our
examples.
Let p = [px , py , pz ]T and v = [vx , vy , vz ]T be the position andvelocity of the centroid of a sub triangle with area A. p andc The Eurographics Association 2010.

1
fdrag = − CD Aeff |vrel |vrel ,2
1
n × vrelflift = − CL Aeff |vrel |(vrel ×),
2|n × vrel |

(15)
(16)

where CD and CL are the drag and lift coefficients, vrel = v −vfluid is the relative velocity of the sub-triangle with respect to∂ηTeff
the fluid. vfluid = [u, u ∂ηx + w z , w] . A is the effective areaof the sub triangle, which depends on the overall structure of
the solid. It is computed as(0
if η(p) < py or n · vrel < 0eff
A =n·vrel( |v | ω + (1 − ω))A otherwise,
rel(17)
where 0 ≤ ω ≤ 1 is a user defined parameter for adjusting
the effective area.
We then add the force fbouyancy + fdrag + flift at position p to thecorresponding solid the sub-triangle belongs to. In the casethe solid is a soft body, we distribute the force to the threevertices of the original triangle weighted by the barycentric
coordinates.
2.3.2. Solids to Fluid
We modify the height and the velocity of the fluid due tosolids using Algorithm 2.
During a time step, a fast moving object may pass throughmany grid cells. In that case the changes to the grid cells andfaces are applied to all the cells along the path (lines 4 to 15).
We also reduce the effect of solids on the fluid exponentiallyNuttapong Chentanez & Matthias Müller / Real-time Simulation of Large Bodies of Waterwith Small Scale Details

Algorithm 2 Solids to Fluid∆t
1: num_substeps = max(1, b|v − vy ŷ| ∆x+ 0.5c)2: Vdisp = n · vrel A∆t
3: sign = (ny > 0)?1 : −1
4: for q = 1 to num_substeps doq∆t5:
ps = p + v

num_substeps(i, j) = closest_gridpoint(p)depth = ηi, j − py
if depth > 0 then9:decay = eλ(−depth)
Vdisp
10:
hi, j += decay num_substepsC(∆x)2 dis
depth
∆t
11:
coeff = min(1.0, decayCadapt ηi, j sign (∆x)
2 A)
12:
ui+ 1 , j += coeff (vx − ui+ 1 , j )2213:vi, j+ 1 += coeff (vz − vi, j+ 1 )2
2
14:
end if
15: end for
6:
7:
8:

with increasing distance to the surface (line 9). The decayrate λ > 0 is set to 1 in all our examples. Cdis > 0 and Cadapt >0 controls how much the solid effects the fluid height andvelocity respectively. We use Cdis = 1.0 and Cadapt = 0.2 inall examples. The height change depends on the amount ofvolume that the triangle sweeps through the fluid during that
Vdisp
sub time-step, num_substeps, divided by the area of the cell,(∆x)2 , see line 10.
2.4. Particle Generation and Simulation
We employ three types of particle in our simulator: spray,splash and foam. Spray and splash particles represent partsof liquid that break free from the height field. Spray particles
are small, fast moving droplets while splash particles represent the remainder. Foam particles represent foam that floats
on the surface of the water. Spray, splash and foam particles
are generated within our simulation framework from breaking waves, waterfalls and interaction of solids with the fluid.
In addition, they can be generated explicitly by user definedsources such as faucets, pouring waters or rain drops. We do
not model bubbles raising from below the water in this work.
Such effects could be added if desired with the method presented in [TSS∗ 07].
2.4.1. Spray, Splash and Foam Simulation
We model spray and splash using particle systems without
particle-particle interaction. The spray and splash particles
need to be generated separately and simulated using different drag coefficients. When a particle is generated froma grid cell, it takes away some mass and momentum fromthe height field. When it falls back, it deposits its mass andmomentum back at another location. If a particle is of raVeff
dius r, the height that will be added/subtract is (∆x)
2 , where

Veff = Cdeposit 43 πr3 is the amplified volume of the particles.
Cdeposit ≥ 1 makes the volume of the particles more thanit physically is. Cdeposit can be used to control the numberof particles active in the simulation. We use the values ofCdeposit between 1 and 10 in our examples. Let ui+ 1 , j be the2nearest velocity along x-axis sample that a particle falls into,
par par par T
parv = [vx , vy , vz ] be its velocity. We can update ui+ 1 , j
2
by using
ui+ 1 , j =ui+ 1 , j hi, j (∆x)2 + vpar
x Veff2

hi, j (∆x)2 +Veff

2

.

(18)

A similar equation is used for wi, j+ 1 . When a splash particle
2
hits the surface we create a foam particle with some probability, depending on the impact speed of the droplet.
We jitter the initial positions of particles by moving them bya random fraction of a time step. This makes them look as
if they were generated somewhere in the middle of the timeinterval. We also add some noise to the initial velocities ofthe particles.
Foam is advected by the velocity field of the fluid simulation and projected onto the fluid surface. Its lifetime is auser-defined parameter modulated with some noise. It nicely
conveys the horizontal swirling water motion in Figure 8.
2.4.2. Breaking Waves
Breaking waves cannot be resolved by a height field model.
In a 2D framework, waves that would break in three dimensions produce disturbing ripples due to numerical instabilities. To solve this problem, we damp such waves and create
particles to make sure we do not lose the effect as shown in
Figure 5a. A cell (i, j) is considered to contain a breaking
wave if it satisfies all of these conditions:1. |∇ηi, j | > αminSplash g∆t∆x// Steep enough to break

prev

2.

hi, j −hi, j∆t> vminSplash // Raising fast and is front of a wave

3. ∇2 ηi, j < lminSplash ,// Is top of a wave

where hprev
i, j is the height in the previous time step. Condition
1 makes sure that the wave is steep enough. It is similar toa condition presented by Thuerey et al. [TMFSG07]. Condition 2 requires that the cell is part of the front of a fast raising
wave. In our examples, we found it to be more robust thanthe condition ∇ηi, j · [ui, j , vi, j ]T < 0 of [TMFSG07]. Condition 3 ensures that we generate particles only near the topof a wave if its wave length is many grid cells wide. We use
vminSplash = 4, αminSplash = 0.45, lminSplash = −4 in our examples. We compute ∇ηi, j using the maximum among the onesided derivatives, to allow detection of a fast height change,as"#∇ηi, j =maxabs (ηi+1, j −ηi, j ,ηi, j −ηi−1, j )∆xmaxabs (ηi, j+1 −ηi, j ,ηi, j −ηi, j−1 )
∆x

,

(19)

c The Eurographics Association 2010.

Nuttapong Chentanez & Matthias Müller / Real-time Simulation of Large Bodies of Waterwith Small Scale Details

Figure 5: a) Breaking wave particles. b) Particles are placed
within the red dotted rectangle.
where maxabs (s,t) = (|s| > |t|)?s : t. We compute ∇2 ηi, j using a central finite differencing as

Figure 6: a) A scene with multiple layers of waterfalls. b)
The waterfall particles are seeded within the dotted red box.
if ηi+1, j < Hi, j thenvwf = [ui+ 1 , j , 0,wi, j+ 1 +wi, j− 12

∇2 ηi, j =ηi+1, j + ηi−1, j + ηi, j+1 + ηi, j−1 − 4ηi, j. (20)(∆x)2

We then generate particles within a rectangle whose edges
are the y-axis and the line perpendicular to the direction of
the gradient passing through the cell center clipped against
cell edges, as shown in Figure 5b. The total volume of theadded particles is proportional to |∇ηi, j |. We jitter these∆xparticle positions by uniform noise sampled form [− ∆x
2 , 2 ].
The x and z components of the particle velocities are givenby the wave velocity [TMFSG07]
and are clamped to the√
maximum water depth,hi, j −h−∇ηi, jprev

g min(hi, j ,havgmax ). The y com|∇ηi, j |ponents is λy ∆t i, j , where we use λy = 0.1 in all examples. We also label a fraction of these particles as spray.

2

2
T

2

]T

pwf
min = [i∆x, Hi, j , j∆x]Tpwf
max = [i∆x + ui+ 1 , j ∆t, ηi, j , ( j + 1)∆x]2else
wi+1, j+ 1 +wi+1, j− 1
2
2 T]
vwf = [ui+ 1 , j , 0,
2
2

T
pwf
min = [(i + 1)∆x + ui+ 1 , j ∆t, Hi+1, j , j∆x]2Tpwf
max = [(i + 1)∆x, ηi+1, j , ( j + 1)∆x]
end ifWe do not subtract the volume of the splash particles fromthe height field because the modified height integration step
explained next already take this into account. The height andvelocity integration is modified as follows:• case 1: Replace ui+ 1 , j with 0 when updating hi+1, j .
2

h

i, j
——–..Replace Eq. 8 with ui+ 1 , j−=g∆t(− ∆x− ∆Hcap ).
2

2.4.3. Waterfalls• case 2: Replace ui+ 1 , j with 0 when updating hi, j .
2

We treat a face (i + 12 , j) as a waterfall face if• case 1:

Hi, j −Hi+1, j> ∆Hcap and Hi, j > ηi+1, j or∆x

• case 2:

Hi+1, j −Hi, j> ∆Hcap and Hi+1, j > ηi, j .
∆x

In each case, the left condition states that the terrain slope islarger than ∆Hcap , where we use 32 in all examples. This indicates an approximate discontinuity in the terrain height. The
right condition makes sure that the lower end of the waterfall is not already filled to the same level as the higher end.
An example with multiple layers of waterfalls is shown in
Figure 6a.
We generate splash particles by sampling uniformly overwf
an axis aligned bounding box defined by pwf
min and pmax aswf
shown in Figure 6b. Their velocities are set to v . Thisbounding box represents the volume of the fluid that would
flow across the face (i + 12 , j) during the current time step.
These variables are computed using the following algorithm:c The Eurographics Association 2010.

hj——–..Replace Eq. 8 with ui+ 1 , j−=g∆t( i+1,
∆x + ∆Hcap ).
2

These modifications have the effect that the cell at the bottom of the waterfall treats the waterfall face as a reflecting
boundary and that the waterfall face’s velocity is updated
specially. The volume that would flows across this boundary is carried by the particles instead. Our waterfall particle
generation algorithm not only enhances the visual richnessof the scene, but also makes the simulation more stable because steep terrain slopes cause numerical instabilities. The
advection of quantities along the velocity field would haveto be modified as well to make sure that they are not traced
up the waterfall. However, we found that this can safely beignored in practice without introducing visual artifacts. Waterfall faces (i, j + 21 ) are treated in a similar manner.
2.4.4. Splashes from Interaction with Bodies
We also generate splashes and spray when rigid and soft
bodies interact with the fluid. For this we go through eachsurface triangle of the solid and check if it sweeps throughNuttapong Chentanez & Matthias Müller / Real-time Simulation of Large Bodies of Waterwith Small Scale Details

the height field during the current time step and its speed isabove a threshold, vthres . If so, we sample the triangle uniformly by recursive subdivision as in Section 2.3. We stop
when the area is roughly the same as the size of splash particles. For each sample, we check whether it sweeps throughthe height field. If so, we generate splash particles. Let ps ,vs , ns be the sample’s position, velocity and normal foundby barycentric interpolation from the original triangle andlet vsfluid be the fluid velocity at ps . The velocity of the particles generated, vpar , is computed using the following algorithm:
vrel = vs − vsfluid , —vrn = (vrel · ns )ns , —vrp = vrel − vrn ,vpar = vsfluid + αrn vrn + αp vrp + αn ns |vrel · ŷ| + αr (vrel − 2(ŷ · vrel )ŷ).

The α’s parameters control the look of the splashes. In our
examples, we use αrn = 0.86, 0.10 ≤ αp ≤ 0.23, 0.00 ≤ αn ≤0.01, and 0.17 ≤ αr ≤ 0.61. We then generate particles at
ps + 2rqns , where 0 ≤ q < num_par, with initial velocity
vpar . We make num_par proportional to |vs | − vthres . The position of these particles are jittered randomly within the radius r.

2.5. Rendering
After a time step, we have a newly updated state of the height
field fluid and the particles. In this section we describe how
to render them. Our goal is to convey the complex fluid flow
features and to enhance the surface details of the fluid.
2.5.1. Wave Crest Alignment and Choppy Waves
The height field fluid can be thought of as a rectangulargrid of quads over the x-z plane. ηi, j determines the ycoordinate of the vertex (i, j) of the grid. Each quad canbe split into two triangles in two ways. Silhouette artifacts are reduced by picking the diagonal that better aligns
with wave crests. We pick diagonal (i, j) − (i + 1, j + 1)
if ηi, j + ηi+1, j+1 > ηi+1, j + ηi, j+1 , otherwise, we pick diagonal (i + 1, j) − (i, j + 1). This test can be implemented
efficiently on the GPU using a geometric shader. We employ this technique in the beach scene (Figure 3) and theopen ocean scene (Figure 2). To improve the appearance of
the water surface further, we make the waves look choppy
as in [LO07], by replacing vertex coordinates (x, z) with(x + c ∂η, z + c ∂η), where c > 0 control the choppiness. We∂x∂z
also clamp the displacement in each direction to not exceed∆x
2 in order to avoid self intersections.
2.5.2. FFT Waves Advection
The height field simulation described so far cannot resolve
waves with wave lengths smaller than the grid resolution ∆x.
With the constraint of being real-time, decreasing ∆x further
is not an option because this would increase the number ofcells and require a smaller time step ∆t for a stable simulation. Instead, we want to add smaller waves on top of theheight field with the following properties:
1. They should be advected with the velocity field.
2. They must not be distorted excessively over time.
3. They disappear if being stretched too much.
4. The method must be relatively cheap.
The key idea is to use an FFT based wave simulation[MWM87,Tes99] with periodic boundary conditions to generate a high resolution wave texture F(s,t) and lookup
this texture for the additional sub-grid displacements of theheight field. The FFT wave texture only contains high frequency waves. The texture coordinates and the weights used
for modulating the additional displacement are computed on
the simulation grid points. They are then bi-linearly interpolated for per-pixel bump mapping or for displacing therendering grid. We use the texture for bump mapping in all
examples. We also use it for displacing a high resolution rendering grid with a grid spacing of 41 ∆x in the boat riding example. Figure 7 shows the comparisons between rendering
without FFT waves, with constant weight FFT waves andwith varying weight FFT waves. Notice how the constantweight case contains unnaturally distorted small waves.
We employ the method present by Max and Becker [MB96]
for the texture coordinates re-generation and texture blending. For each grid cell (i, j) we additionally store 3 setsof texture coordinates (ski, j ,ti,k j ), k = 1..3 which are initialized with (i∆x, j∆x). They are transported with the velocity
field using semi-Lagrangian advection. We reset (ski, j ,ti,k j )to (i∆x, j∆x) every τ seconds with a phase shift of τ/3 between k and (k + 1) mod 3. The weight used for multiplying the displacement ωi, j should depend on how much thevelocity field would stretch the wave. We use ωi, j = e−Ωµi, j ,where µi, j = max3k=1 µki, j . µki, j is the absolute value of theeigenvalue with maximum magnitude of the Green strainof (sk ,t k ) at (i, j). The Greenstrain can be computed asT12 (DD − I), where D =

∂sk
∂x
∂t k
∂x

∂sk
∂z
∂t k
∂z

. Ω > 0 controls how

fast the small waves disappear with stretching, for which
we use values between 0.5 to 2 in our examples. The final displacement for a given position is ω ∑3k=1 wk F(sk ,t k ),where ω, sk ,t k are bi-linearly interpolated from the grid,t−t kwk = 13 (1−cos(2π τ 0 )), t0k is the last time the set k is regenerated. The phase shift ensures a constant blending weight.
The choice of τ trades off the distortion of the texture coordinates against the coherency of the waves. We use τ ≤ 3seconds in our examples. One could also adapt a texture coordinate regeneration method presented by Neyret [Ney03]
to not have to make this tradeoff, but it is more expensive.
c The Eurographics Association 2010.

Nuttapong Chentanez & Matthias Müller / Real-time Simulation of Large Bodies of Waterwith Small Scale Details

Figure 7: Left: No FFT Mid: FFT constant weight (Ω = 0) Right: FFT varying weight (Ω = 1)
2.5.3. Particles Rendering
We render the splash particles using the screen space technique presented by Van der Laan et al. [VdLGS09] with themodification of using a bi-lateral filter with support radiusdefined in world space [RD10], so that the fluid appearance
is view-independent. For the final rendering of the particle
surface, we use the same pixel shader as the one for rendering the height field.
Spray particles are rendered as an elongated ellipse along
the direction of their velocity to emulate the motion blur effect as in [VdLGS09]. Foam particles are rendered as diffuse
disks with normals perpendicular to the height field water
surface.

We are currently experimenting with a method to simulate
this phenomena. Another area of future work is in renderingof the water represented by the height field and the particles
seamlessly and efficiently.

Grid
Par

Boat
WaterF PML
Beach Ocean
900x135-128x128 -128x128 -128x128-256x256
250K
56K
2K
220K
83K

Table 1: Grid size and average number of foam, spray, andsplash particles active in a given frame.

3. Results and Discussion
We implement our method both on the CPU as well as on the
GPU using CUDA. All timings are done on an Intel Core i7,
2.67 GHz with 4GB of memory and NVIDIA GTX480. We
use NVIDIA’s PhysX SDK for rigid body, soft body, andcloth simulation. The size of the grid and the number of particles used for each example are stated in Table 1. We use∆t = 16.66ms in all examples, except the boat riding example, where we use ∆t = 20ms. Table 2 show the timings forvarious examples. For scenes indicated with the letter G, all
the steps are executed on the GPU except for the solid coupling. Porting the solid coupling part to the GPU is part of
future work. For the CPU cases, the height field simulationand the particles simulation are multi threaded.
A drawback of our method is that although the height integration step conserves volume, the coupling and the overshooting reduction do not. However, we found that volume
loss is less crucial in these situations and a global mass adjustment technique such as the one proposed by [Hes07]
could be employed if desired. Another drawback is that our
stability enhancements cannot guarantee unconditional stability. However, they allow us to explicitly integrate with areasonable time step in all of our examples.
When a dense object falls into the water and fully submerges,
there is usually a secondary splash called Rayleigh jet due tothe collapse of water over the cavity formed by the object.
c The Eurographics Association 2010.

BoatC1
BoatC2
BoatC3
BoatC4
BoatG
WaterFG
PMLG
BeachG
OceanG

Total
94.50
76.75
73.88
69.75
18.05
7.97
3.78
9.88
13.13

HF
42.25
25.69
23.25
19.06
2.19
0.75
0.78
0.75
1.28

Gen
5.11
5.11
5.11
5.09
1.36
0.82
0.13
1.87
1.16

Par
16.18
9.91
8.18
7.55
1.15
0.44
0.19
0.90
0.56

Cou
2.46
2.31
2.24
2.31
2.32
0.00
0.00
0.00
3.13

Table 2: Timing for various examples in milliseconds. Total is the frame time including rendering. HF includes advection, height integration, velocity integration, PML, overshooting reduction and FFT texture coordinate advection.
Gen includes all particle generation steps. Par includes particles simulation. Cou is the coupling of solids with fluid. Thelines BoatC1...C4 contain the timings of the boat demo running on the CPU with 1-4 threads. The other lines show timings on the GPU. PML shows the all absorbing boundaries
case. Beach shows the timings for the case with overshooting
reduction.

Nuttapong Chentanez & Matthias Müller / Real-time Simulation of Large Bodies of Waterwith Small Scale Details

References[ATBG08] A NGST R., T HÜREY N., B OTSCH M., G ROSS M.:
Robust and Efficient Wave Simulations on Deforming Meshes.
Computer Graphics Forum 27 (7) (October 2008), 6, 1895 –1900. 2, 3
[BAB09] B URRELL T., A RNOLD D., B ROOKS S.: Advectedriver textures. Comput. Animat. Virtual Worlds 20, 2–3 (2009),
163–173. 2[Pea86] P EACHEY D. R.: Modeling waves and surf. SIGGRAPH
Comput. Graph. 20, 4 (1986), 65–74. 2[RD10] ROUSLAN D IMITROV M. S.: SPH fluid rendering, in
submission to siggraph talk, 2010. 9[SFK∗ 08] S ELLE A., F EDKIW R., K IM B., L IU Y., ROSSIGNAC
J.: An unconditionally stable MacCormack method. J. Sci. Comput. 35, 2-3 (2008), 350–371. 3[Bri05] B RIDSON R.: Shallow water discretization, Lecture notes
Animation Physics. University of British Columbia, 2005. 3[SM09] S ÖDERSTRÖM A., M USETH K.: Non-reflective boundary conditions for incompressible free surface fluids. SIGGRAPH
Talks (2009). 4[Bri08] B RIDSON R.: Fluid Simulation for Computer Graphics.
A K Peters, 2008. 2[SP09] S OLENTHALER B., PAJAROLA R.: Predictive-corrective
incompressible sph. In Proc. SIGGRAPH (2009), pp. 1–6. 2[CL95] C HEN J. X., L OBO N. D . V.: Toward interactive-rate
simulation of fluids with moving obstacles using navier-stokes
equations. Graph. Models Image Process. 57, 2 (1995), 107–116.
2[TDG00] T HON S., D ISCHLER J.-M., G HAZANFARPOUR D.:
Ocean waves synthesis using a spectrum-based turbulence function. In Proc. CGI (2000), p. 65. 2[CS09] C ORDS H., S TAADT O.: Real-time open water environments with interacting objects. In Proceedings of Eurographics
Workshop on Natural Phenomena (EGWNP) (2009). 2[FR86] F OURNIER A., R EEVES W. T.: A simple model of ocean
waves. SIGGRAPH Comput. Graph. 20, 4 (1986), 75–84. 2[Hes07] H ESS P.: Extended Boundary Conditions for Shallow
Water Simulations. Master’s thesis, ETH Zurich, 2007. 2, 3,
9[HHL∗ 05] H AGEN T. R., H JELMERVIK J. M., L IE K.-A.,
NATVIG J. R., H ENRIKSEN M. O.: Visual simulation ofshallow-water waves. Simulation Modelling Practice and Theory 13, 8 (2005), 716–726. 2[Tes99] T ESSENDORF J.: Simulating ocean water. SIGGRAPH
course notes, 1999. 2, 8
[TMFSG07] T HUREY N., M ULLER -F ISCHER M., S CHIRM S.,
G ROSS M.: Real-time breakingwaves for shallow water simulations. In Proc. Pacific Conf. on CG and App. (2007), pp. 39–46.
2, 6, 7
[TSS∗ 07] T HÜREY N., S ADLO F., S CHIRM S., M ÜLLER F ISCHER M., G ROSS M.: Real-time simulations of bubbles andfoam within a shallow water framework. In Proc. SCA (2007),
pp. 191–198. 6[vBBK08] Š T́AVA O., B ENEŠ B., B RISBIN M., K ŘIVÁNEK J.:
Interactive terrain modeling using hydraulic erosion. In Proc.
SCA (2008), pp. 201–210. 2[HNC02] H INSINGER D., N EYRET F., C ANI M.-P.: Interactive
animation of ocean waves. In Proc. SCA (2002), pp. 161–166. 2[VdLGS09] VAN DER L AAN W. J., G REEN S., S AINZ M.:
Screen space fluid rendering with curvature flow. In Proc. I3D(2009), pp. 91–98. 9[HW04] H OLMBERG N., W ÜNSCHE B. C.: Efficient modelingand rendering of turbulent water over natural terrain. In Proc.
GRAPHITE (2004), pp. 15–22. 2[WMT07] WANG H., M ILLER G., T URK G.: Solving general shallow wave equations on surfaces. In Proc. SCA (2007),
pp. 229–238. 2, 3

[Joh08] J OHNSON S. G.: Notes on perfectly matched layers. online mit course notes, July 2008. 4[YHK07] Y UKSEL C., H OUSE D. H., K EYSER J.: Wave particles. In Proc. SIGGRAPH (2007), p. 99. 2, 5

[KM90] K ASS M., M ILLER G.: Rapid, stable fluid dynamics forcomputer graphics. In Proc. SIGGRAPH (1990), pp. 49–57. 2[YNBH09] Y U Q., N EYRET F., B RUNETON E., H OLZSCHUCH
N.: Scalable real-time animation of rivers. Proc. Eurographics
28, 2 (2009). 2[LO07] L EE R., O’S ULLIVAN C.: A fast and compact solver
for the shallow water equations. Virtual Reality Interations andPhysical Simulation 1 (2007), 51–58. 2, 8
[LvdP02] L AYTON A. T., VAN DE PANNE M.: A numericallyefficient and stable algorithm for animating water waves. The
Visual Computer 18, 1 (2002), 41–53. 2, 3
[MB96] M AX N., B ECKER B.: Flow visualization using moving
textures. In Proc. LCAS/LaRC Symposium on Visualizing TimeVarying Data (1996), pp. 77–87. 8[MFC06] M AES M. M., F UJIMOTO T., C HIBA N.: Efficient animation of water flow on irregular terrains. In Proc. GRAPHITE(2006), pp. 107–115. 2[MWM87] M ASTIN G. A., WATTERBERG P. A., M AREDA J. F.:
Fourier synthesis of ocean scenes. IEEE Comput. Graph. Appl.
7, 3 (1987), 16–23. 2, 8
[MY97] M OULD D., YANG Y.-H.: Modeling water for computer
graphics. Computers & Graphics 21, 6 (1997), 801–814. 2

Appendix
3.1. PMLs for z-moving waves
hi, j += (−γi, j (hi, j − hrest ) + ψi, j )4t
1
wi, j+ 1 += − (γi, j+1 + γi, j )wi, j+ 1 4t
2
2
2
(ui+ 1 , j − ui− 1 , j )22
ψi, j += −λupdate ψi, j
4t
4x
ψi, j ×= λdecay

(21)
(22)
(23)
(24)

We evolve the field ψ only in the regions where we wantto damp a z-moving waves. γ plays the same role as σ in
Section 2.1.4.

[Ney03] N EYRET F.: Advected textures. In Pro. SCA (2003),
pp. 147–153. 8[OH95] O’B RIEN J. F., H ODGINS J. K.: Dynamic simulation ofsplashing fluids. In Proc. Computer Animation (1995), p. 198. 2
c The Eurographics Association 2010.

Nuttapong Chentanez & Matthias Müller / Real-time Simulation of Large Bodies of Waterwith Small Scale Details

Figure 8: A boat driving on an infinitely long river that flows over an irregular terrain. This scene demonstrates all techniquespresented in the paper: Small waterfalls, swirling vortices, foam, spray, splashes, small scale waves, interaction of water withrigid and soft bodies. The simulation grid is shifted in whole multiple of grid spacing ∆x to be approximately centered aroundthe boat at all time.

c The Eurographics Association 2010.