Fluid Dynamics with Erosion


I really wanted to do something with fluid dynamics for my final project. I wanted to be able to pour a pitcher of water on a bank of packed sand and watch the formation of eroded channels and alluvial fans where sand is deposited. This is an interesting problem because the fluid flow affects its boundaries through erosion and deposition, which in turn affect the fluid flow. To my knowledge there had been no work in computer graphics on small-scale erosion phenomena of this nature. There has been work done on erosion of valleys and river beds over large time scales. Most of these approaches use empirical models for sediment transport and fluid flow. In order to get visually interesting results I wanted to solve for fluid flow with the Navier-Stokes equations for incompressible fluids and combine that with a physically-based sediment transport model.

There has been much work done in simulating sediment transport for engineering applications. Numerical simulation is often performed before building dams, weirs, bridges and other structures  to predict and avoid possible sedimentation and scouring around the structures. The focus of the models used for these simulation is physical accuracy rather than visual effect. These models form the basis for the model I used in this project.

Fluid Model

The Navier-Stokes equations provide an accurate model of laminar fluid flow. For more details on these equations and their simulation in computer graphics, look here.

The fluid solver I am using for this project is based on [Griebel et al. 1998]. The Navier-Stokes equations are discretized on a staggered grid. After computing a tentative velocity field, a pressure correction term to make the field divergence free is computed using successive over-relaxation (SOR). Boundaries are assumed to lie along the grid edges at boundary cells. I never got around to trying the fluid solver used by [Stam 1999]. All of the necessary components are already  in the code because I use that model for suspended sediment transport, but I ran out of time. In particular, I would like to compare the amount of numerical dissipation that occurs from the semi-lagrangian advection step used in Stam's solver.

Sediment Transport Model

There are two modes of sediment transport. Bed-load transport occurs when particles of sediment slide, roll, or saltate along the bottom. These particles move along in a region above the bed called the bed-load layer. The bed-load layer has a height corresponding to the maximum height achieved by a saltating particle. When the velocity of the fluid is high enough particles begin to move into the second mode of transport, suspended transport. Suspended particles entrained at the top of the bed-load layer move through the fluid by advection and diffusion.

There are a large number of empirical models for sediment transport (see [Chanson 1998]) which do not necessarily agree. Since I was more concerned about getting the visual aspect of sediment transport than physically accurate simulation I decided to start with a simplified model. I did away with the complicated expressions for coefficients based on physical parameters and replaced them with a single coefficient that could be set with a slider to achieve the desired effect. The model does not include bed-load transport. Though the contribution of bed-load transport to the total sediment transport is not neglible, suspended transport is the most important factor.

The main idea behind the model I used is that the capacity for a fluid to carry suspended sediment is related to the velocity of the fluid and the fall-rate of the sediment. A sediment particle has a terminal fall velocity that depends on its shape, size, and density relative to the fluid. Smaller particles such as clay and silt tend to have smaller fall velocities and thus remain suspended longer, whereas gravel and sand have much higher fall velocities. The critical velocity at which sediment begins to move into suspension is related to the fall-velocity. A smaller disturbance in the fluid is required to suspend particles with a low fall velocity. Similarly, the sediment capacity is also related to fall velocity. When the sediment capacity above the bed is higher than the concentration of sediment there is a net flux from the bed into the fluid - erosion. When the concentration is higher than the sediment capacity, some of the sediment is deposited. I add a term that increases the sediment capacity and critical velocity according to the slope of the terrain.

Once the sediment has entered into suspension it is transported using the techniques described in [Stam 1999]. The semi-Lagrangian advection uses the velocity field to take an Eulerian step backward in time to determine what the concentration would have been advected to the current position. In general, this step does not conserve mass, especially in the presence of vorticies. A more serious problem occurs at the bottom boundaries. The sediment moves downward according to its fall velocity. On the row of cells above the bottom boundary, the advection step looks backward and does not see the wall. All sediment that was in the bottom row is lost because it is not advected anywhere. To deal with this problem I advected the sediment into the top row of boundary cells and added the advected sediment back into the row above the boundary. This almost conserves volume. Volume actually goes up slightly with time. It has the disadvantage that the time-step has to be limited such that the sediment does not move more than one cell downward in a single step. Since the fall velocity is usually very small, this hasn't been a significant limitation.

As the sediment is eroded and deposited, the terrain changes. I represent the terrain as a heightfield. For the purposes of the fluid simulation, the heightfield is approximated with boundary cells. Velocity values at the boundary cells are used to satisfy boundary conditions. No boundary cells are allowed that have fluid on opposite sides as the value stored at the boundary cell would have to be shared by both edges. At each iteration I update the heightfield according to the sediment transport. Then I build a set of boundary cells that conform to the new heightfield. I then delete boundary cells that have fluid on opposite sides. The velocities and pressures in the cells are not updated until the final status of the cell has been determined.


I did not get as far as I had hoped on this project. I was never able to get a free surface for the fluid. The preliminary results are very nice. The erosion and deposition behaves plausibly. The evolving bed forms affect the fluid flow.

I ran the simulation on a 150x40 grid with periodic boundaries to simulate a river bed. The user can apply velocity and add sediment concentration with the mouse. The code is completely unoptimized yet the simulation easily runs at 25 fps on a 1.7 GHz Pentium 4 with a GeForce 3 graphics card.  Here are some images produced by the simulation.

Plumes of sediment and a crater left after a blast of fluid towards the bed.

Erosion off the tops of the dunes in a channel flowing to the right.

Grid and boundary cells used for fluid computation.

Strange spiky artifact that sometimes occur.

These spikes are a result of a defect in my model. I think what is happening is that sediment moves out of suspension near the tops of the spikes leaving little left by the time it falls to the bottoms. Over time the spikes grow taller and taller, increasingly depleting the sediment reaching the troughs. One possible solution to this would be to add instability based on slope. Slopes higher than a given threshold become unstable and redistribute their sediment to neighbors.

Future Work

The results of this project are extremely encouraging. There are many things that I would like to do to improve it. First of all, I would like to implement a more sophisticated transport model involving bed-load transport. I would like to experiment with different methods for suspended transport that conserve mass. I would like to rigorously test and compare different fluid solvers. The solver could be improved to respond better to the actually shape of the heightfield. This would help eliminate blocky artifacts that are sometimes visible due to the boundary cell approximation. Adding a free surface would make the simulation much more interesting as bed-forms would affect the formation of surface waves. Heterogeneous sediments would allow for the formation of layers that are often seen in sedimentary rocks. Finally I would like to implement portions of the simulation on the GPU to facilitate larger grid sizes. I plan to eventually make this into a screen-saver type thing that I will release to the public.

Related Material

Proposal (Powerpoint)
Progress Report  - 4/07/03
CFD Lecture (Powerpoint)
Sediment Transport / Final Presentation (Powerpoint)

Video 54s (DivX - 5.3 MB)


Sediment Transport

Chanson, H. 1999. The Hydraulics of Open Channel Flow: An Introduction. Arnold.

Chiba, N., Muraoka, K., and Fujita, K. 1998. An erosion model based on velocity fields for the visual simulation of mountain scenery. The Journal of Visualization and Computer Animation, vol. 9, 185-194.

Haupt, B. J., Seidov, D. and Stattegger, K. 1999. SEDLOB and PATLOB: Two numerical tools for modeling climatically forced sediment and water volume transport in large ocean basins. In Computerized Modeling of Sedimentary Systems. Springer-Verlag, Berlin.

Kelley, A., Malin, M., and Nielson, G., 1988. Terrain simulation using a model of stream erosion. In Computer Graphics (Proceedings of SIGGRAPH 88). vol. 22, 263-268.

Musgrave, K., Kolb, C., and Mace, R., 1989. The synthesis and rendering of eroded fractal terrains. In Computer Graphics (Proceedings of SIGGRAPH 89). vol. 23, 41-50.

Nagashima, K. 1998. Computer generation of eroded valley and mountain terrains. The Visual Computer, vol. 13, 456-464.

Wu, W., Rodi, W. and Thomas, W. 2000. 3D numerial modeling of flow and sediment transport in open channels. Journal of Hydraulic Engineering, 4-15.


Chen, J., and Lobo, N. 1994. Toward interactive-rate simulation of fluids with moving obstacles using the navier-stokes equations. Computer Graphics and Image Processing, 107–116.

Chen, S., Johnson, D., Raad, P. and Fadda, D. 1997. The surface marker and micro cell method. International Journal of Numerical Methods in Fluids, 25, 749-778.

Enright, D., Marschner, S. and Fedkiw, R. Animation and rendering of complex water surfaces,   In Proceedings of ACM SIGGRAPH, 736-743.

Griebel, M., Dornseifer, T., and Neunhoeffer, T. 1998. Numerical Simulation in Fluid Dynamics: A Practical Introduction. SIAM Monographs on Mathematical Modeling and Computation. SIAM.

Foster, N., and Metaxas, D. 1996. Realistic animation of liquids. Graphical Models and Image Processing, 471–483.

Foster, N., and Fedkiw, R. 2001. Practical animation of liquids. In Proceedings of SIGGRAPH 2001, 23–30.

Kass, M., and Miller, G. 1990. Rapid, stable fluid dynamics for computer graphics. In Computer Graphics (Proceedings of SIGGRAPH 90), vol. 24, 49–57.

O’Brien, J., and Hodgins, J. 1995. Dynamic simulation of splashing fluids. In Proceedings of Computer Animation 95, 198–205.

Stam, J. Stable fluids. 1999. In Proceedings of SIGGRAPH 99, 121-128.

Tome, M. F.  and McKee, S.  1994. GENSMAC: A computational marker-and-cell method for free surface flows in general
domains. Journal of Computational Physics, 171-186.