This page does not represent the most current semester of this course; it is present merely as an archive.

Hydraulic erosion is a dominant force in shaping most landscapes on Earth. The basic process runs as follows:

  1. Rain falls from the sky, putting water on the ground
  2. Water flows down hill
  3. Water can store and transport sediment
    • Faster water can store more than slower water
    • Faster water erodes more, removing material from the ground
  4. Water both evaporates and is absorbed by the ground

Within these general guidelines there are many additional factors to consider. Water has momentum and interacts with other water; water dynamics change with temperature and amount of carried sediment; the ground is made of many types of material, some more easily eroded than others; drying out can change the material properties of the ground; and so on.

Simple hydraulic erosion works as follows:

  1. An initial non-eroded terrain is created, either by a fractal method like faulting planes or by an artist-supplied base mesh.

  2. A few thousand iterations of water moving across the terrain are simulated. There are two broad ways to do this:

    1. Water is modeled as a grid of water heights and other properties. Basic laws of water flow are used to update this grid. A few thousand iterations are needed to get reasonable erosion effects, each of which processes the entire grid.

    2. Water is modeled as a sequence of individual particles. Particles are distributed based on rainfall models and roll down hill. Hundreds of thousand of iterations are needed to get reasonable erosion effects, each of which processes the path of just one particle.

This page explores a simple version of each of these methods.

1 Grid-based approach

A simple grid-based erosion technique was published by Musgrave, Kolb, and Mace in 1989. It runs as follows:

  • Let each vertex have an altitude aa, water volume ww, and suspended sediment ss.

    The top of the water at a vertex is thus at height w+aw+a; that’s important because the difference in top-of-water heights is what decides which way water flows (and how quickly).

  • Move water and sediment between adjacent vertices ii and jj using Δw=min(wi,(wi+ai)(wj+aj))\Delta w = \min\big(w_i, (w_i+a_i)-(w_j+a_j)\big)

  • If Δw0\Delta w \le 0, increase aia_i and decrease sis_i by KdsiK_d s_i, where KdK_d is a deposition rate constant you pick

  • Otherwise,

    • move Δw\Delta w from wiw_i to wjw_j
    • let c=KcΔwc = K_c \Delta w, where KcK_c is a sediment carrying capacity constant you pick
    • If si>cs_i > c
      • add cc to sjs_j
      • add Kd(sic)K_d(s_i-c) to aia_i
      • set sis_i to (1Kd)(sic)(1-K_d)(s_i-c)
    • otherwise
      • add si+Ks(csi)s_i+K_s(c-s_i) to sjs_j, where KsK_s is a soil softness constant you pick
      • subtract Ks(csi)K_s(c-s_i) from aia_i
      • set sis_i to 00

In the original paper this algorithm was said to apply to each neighboring vertex with this caveat:

In a full two-dimensional implementation, one must take care to distribute water and sediment to all neighboring lower vertices in amounts proportional to their respective differences in overall elevation.

This not only requires some care to pick the distributions, but also requires care so that the order in which you compute the motion from cell ii to jj vs the motion from cell jj to kk does not change the results.

Various realizations of this distribution criteria has been used; one of the simplest is a two-pass system

  1. In the first pass, find how water and sediment should move:
    • Pick Δw\Delta w for each pair of cells such that the total Δw\Delta w out of cell ii is at most wiw_i and is distributed between all neighbors with aj+wj<ai+wia_j+w_j < a_i+w_i proportionally to how different those are
    • Accumulate the corresponding changes to ww, aa, and ss but do not apply them yet
  2. In the second pass, apply the changes to ww, aa, and ss

Parameters in this method:

  • how much ww to add initially (rainfall amounts)
  • how often to have more ww added (rainfall frequency)
  • whether to have ww reduce a bit each step (evaporation) or not, and if so by how much
  • KdK_d, KcK_c, and KsK_s controlling erosion rates

2 Particle-based approach

The first particle-based erosion technique was published by Chiba, Muraoka, and Jujita in 1999. That paper described a fairly involved algorithm, including a model of the under-cutting effects of turning rivers, though their simulations did not show that complexity’s impact.

A simpler version works as follows:

  1. Pick a random point and drop a particle there with
    • Zero velocity v\vec v
    • Some initial water volume ww (i.e. a positive starting value you pick)
    • No sediment ss
  2. While the particle retains volume,
    • Accelerate the particle by adding KaK_a times the non-vertical component of the surface normal to its velocity
      • KaK_a is an acceleration constant you pick
      • you’ll need to dynamically compute the normal based on the ever-changing heights (and normalized it after computation)
    • Slow the particle by multiplying its velocity by (1Kf)(1-K_f)
      • KfK_f is a friction constant you pick
    • Move the particle by adding its velocity to its position
      • if it goes off the map, stop working with this particle
    • Compare the sediment it is carrying (ss) to the sediment it could carry (KcvwK_c \|\vec v\| w)
      • KcK_c is a sediment carrying constant you pick
      • if it has more sediment than it should, move KdK_d of the difference from ss to the altitude of the terrain particle
        • 0Kd10 \le K_d \le 1 is a deposition constant you pick
      • if it has less sediment than it should, move KsK_s of the difference from the altitude of the terrain to ss
        • 0Ks10 \le K_s \le 1 is an soil softness constant you pick
    • reduce ww by some small fixed evaporation rate you pick
      • if w0w \le 0, stop working with this particle
  3. Repeat the above several hundred thousand times

The above has a particle moving in analog steps across a discrete grid. While it is somewhat more accurate to interpolate particle positions onto the nearby terrain vertices, it is generally adequate to simply round to the nearest vertex instead when finding normals and lifting and depositing sediment.

Parameters in this method:

  • initial ww per particle and evaporation rate (together giving a lifetime of the particle)
  • KaK_a and KfK_f controlling particle motion
  • KdK_d, KcK_c, and KsK_s controlling erosion rates