Around The World, Part 23: Hydraulic erosion

Fri, 6 Jun 2025 by Thomas ten Cate · Comments
 Game Development, Around the World

As I mentioned last time, I’m currently working on a full rewrite of the game, with a focus on building a solid technical foundation first. But because much of that is boring work, I allowed myself a fun side quest: hydraulic erosion.

In nature, erosion is a hugely important factor in shaping landscapes. And arguably the most significant form of erosion is hydraulic erosion, that is, erosion by flowing water. It results in typical branching folds in the landscape:

Credit: OutcropWizard via Wikimedia Commons, CC-BY-SA 4.0

Just so we’re on the same page, I’ll first describe my understanding of the hydraulic erosion process. Rain falls on the ground and flows downhill. It combines to form streams, which merge to form rivers, which merge into bigger rivers, until the water eventually flows into the sea. Along the way, it wears down the terrain that it flows over, lowering it. These loose rocks, pebbles and grains of sand, collectively known as sediment are carried downstream. Where the terrain flattens out, the water flows less quickly, and the carried sediment settles down on the bottom, increasing the height of the terrain there.

It’s difficult to recreate the effects of hydraulic erosion directly through equations, although this amazing Shadertoy by Felix Westin does a very impressive job, at the expense of being pretty complicated. Instead, most people resort to some kind of simulation of the erosion process over multiple time steps, simulating water flow, sediment pickup and deposition, and evaporation. Those simulations can get pretty complex and time consuming, so I thought: maybe there’s a simpler way?

Attempt 1: blurosion

It’s always worth starting out with “the simplest thing that could possibly work”. Because it’s simple, it doesn’t take much time to implement, and it may be good enough. Even if it isn’t, it’ll teach you something valuable.

With that in mind, I tried the following idea. The steeper the terrain is, the faster the water will flow, and thus the more sediment it will carry downhill. So here’s the simplest possible erosion algorithm I could think of: for each grid cell, find its lowest neighbour, and compute the difference in height. Reduce my own height by some fixed fraction of that difference (erosion), and increase the neighbour’s height by the same amount (deposition). Repeat this a couple of times.

For example, suppose the current cell’s height is 12, its lowest neighbour has height 4, and we’ve configured the erosion strength at 25%. The height difference is 8, so we’ll want to move 25% × 8 = 2 units of material. The current cell’s height then becomes 12 - 2 = 10, and the neighbour’s height is increased to 4 + 2 = 6.

Let’s try it on this height map, which is mostly just simplex noise with some tectonic boundary effects on top:

And here’s the result after some iterations of the above algorithm:

It looks like I just invented a rather slow blur algorithm instead. Not quite what I’d hoped for. In fact, it’s very similar to the thermal erosion algorithm I described over two years ago and apparently forgot about!

But why didn’t it produce these typical ridges and valleys? Here comes the lesson: hydraulic erosion is, in some sense, self-reinforcing. If by random chance a slight valley forms, it’ll collect slightly more water from its surroundings, causing it to erode more quickly, and becoming more pronounced, which will cause it to collect even more water, and so on. In technical terms, the amount of erosion that happens to a point is strongly dependent on the drainage area of that point, which is the total area that’s upstream from that point. My simplistic idea did not capture that.

Attempt 2: simulation

Looking at the literature, I quickly found a well-known paper by Xing Mei et al., Fast Hydraulic Erosion Simulation and Visualization on GPU. I’m not running my algorithms on the GPU at the moment, but that doesn’t matter; the algorithm works equally well on a CPU and can easily be parallelized.

Mei models water flow as a column of water on top of each terrain cell, connected to its neighbours by virtual pipes. Water flow through these pipes is simulated, which results in a flow velocity field. The sediment transport capacity is then calculated for each cell, depending on the slope and velocity. Sediment is then picked up or deposited to bring the amount of dissolved sediment towards this target capacity. Finally, sediment is moved (advected) through the flow velocity field.

The images in the paper aren’t very impressive, which may in part be due to limitations of the algorithm. Several people have taken the ideas and built upon them, for example Dylan Mcleod et al. in Interactive Hydraulic Erosion Simulator and Lan Lou in Terrain erosion sandbox in WebGL. Those people do get good results, so I knew it was possible.

However, I failed to get good output from Mei’s algorithm. I’m not sure why – maybe my implementation had bugs, maybe the combination of (about ten) parameters was wrong, but I always got weird artifacts at best, and numerical instability at worst. Here’s the best result I managed to get:

It does show branching patterns, but they have a strong tendency to be strictly horizontal or vertical, which gives a very artificial look. You may need to zoom in to see it.

There may have been ways to fix this, but the algorithm was already complex enough. So I started shopping around for something less complicated, but still sophisticated enough to do what I wanted.

Attempt 3: stream power law

While stumbling around, I found the blog post Simulating worlds on the GPU by David A Roberts. In it, he mentions the stream power law, which is a very simple empirical equation that models how quickly terrain erodes:

E = K * A^m * S^n

where

• E is the rate of erosion,
• A is the drainage area,
• S is the slope,
• K, m and n are constants.

A simple equation just two inputs, drainage area and slope, and only three parameters to tune (maybe just two because n should be about twice as large as m) – but it still involves that all-important factor, the drainage area. Did I hit gold?

Slope is easy enough to compute; it’s just the length of the gradient vector. However, how do we find the drainage area? This is where I found another paper, by Schott et al. (2002), Large-scale terrain authoring through interactive erosion simulation. They do a lot of work that I’m not interested in here, but also happen to present an algorithm to approximate the drainage area iteratively. It’s rather simple:

• Set the drainage area of each cell to 1.
• Repeat for some number of iterations:
  • Set the drainage area of each cell to 1 (the cell itself) plus the drainage area of all neighbouring cells that drain into the current cell (i.e. that have the current cell as their lowest neighbour).

This will eventually converge to the true drainage area, but we can stop it sooner if we want.

A complication is that the drainage area depends on the shape of the terrain, and erosion will change that shape. So I decided to iteratively erode the terrain based on the stream power law in the same loop that also updates the drainage area.

Does it work?

It still has directional (horizontal/vertical) artifacts, but it’s getting somewhat acceptable.

However…

As you can see, the algorithm produces a large number of inland lakes below sea level. There are also many higher-up regions that are completely surrounded by higher terrain. These are known as endorheic basins, and are quite rare in nature. The way I understand it, there are two reasons for this. Firstly, they would fill up with water, which carries sediment, eventually raising the basin floor until the point where it is no longer a basin. And secondly, because the water eventually overflows the edge of the basin, it will start carving out a path for the basin to drain into, turning it into a regular lake with an outflow, and eventually maybe just a riverbed.

It makes sense that the algorithm produces basins, because they have a large drainage area, causing their bottom to sink. What the algorithm fails to capture, though, is that the stream power law only applies to flowing water, and the water in a basin is standing still.

I considered various approaches of detecting basins, and then either filling them or carving a way out, but decided that it would be too fiddly to get right, and possibly also too slow. Time for another tack.

Interlude: a better starting point

I decided to turn off the hydraulic erosion algorithm for the moment, and apply some long-wanted improvements to what I already had. This is basically just a rehash of previous stuff that I implemented in C# on a sphere, but now done in Rust on a plane.

Firstly, we know that water runs downhill and eventually reaches the sea. So, very broadly speaking, the closer to the coast a point is, the more water it receives, the more it erodes, and the lower it tends to be. So I computed the shortest distance to the coast at every point, mapped this through a configurable curve, and applied that to the terrain as a starting point:

Of course, for various reasons, the coast isn’t equally steep everywhere. So let’s use some low-frequency simplex noise to modulate the distance by, before mapping it to a height:

That’s better. Now we have steep coasts, as well as flat regions that reach farther inland (hello Netherlands!).

Next up was to improve the noise. What I had previously was just some octaves of simplex noise, which I’d hoped to make more interesting through the erosion algorithm:

(On top of the simplex noise, you’ll also notice some mountain ranges added by tectonic boundaries, which I blogged about extensively before.)

Inspired by the article on value noise derivatives by Inigo Quilez, I modified the noise according to a simple idea: mountains have steeper slopes, and also more rugged terrain. In terms of simplex noise octaves, this translates to: if the gradient of the noise is larger, subsequent smaller octaves of noise get bigger amplitude. I implemented a configurable curve so I can control this effect, which I call “slope dependent amplitude”:

With this curve applied, we can see that flatter regions become smoother, and steeper regions become rougher:

It makes the terrain much more varied and interesting.

Because the slope is apparently not continuous, it also results in some weird artifacts in the form of straight lines… but I have an old trick up my sleeve: domain warping. If we modify the input point of the noise by yet another layer of noise (or rather two, one for x and one for y), those straight lines get broken up nicely, and the overall terrain becomes even more interesting:

This could almost pass for something produced by erosion, without needing an actual erosion algorithm, but I wasn’t satisfied yet. There was one more idea I wanted to try.

Attempt 4: droplets

So far, I’ve only talked about grid-based erosion algorithms. But there is another category, one that is at least as popular: particles. The idea is to simulate rainfall one “droplet” at a time, tracing its path downhill, picking up and depositing sediment as it moves over the terrain. Do that for enough random droplets, and you get quite a convincing landscape. For a more thorough introduction, check out the amazing video Coding Adventure: Hydraulic Erosion by Sebastian Lague. It’s so good I’ll just embed it here:

However, there is one major drawback to this approach: it’s relatively slow. To get some decent results, you need to simulate many thousands of droplets, each for several dozen time steps. Because each droplet affects the terrain in random places, it’s hard to parallelize this across multiple cores, and it’s also not very friendly to CPU caches.

Another drawback is that a droplet can affect the terrain far away from where it’s initially spawned. This makes it difficult to apply erosion to individual terrain chunks, something I’ll want to do later.

For these reasons, I initially didn’t consider using such a droplet algorithm. But maybe I could take some ideas from it, and from the grid-based methods, and come up with some hybrid approach that combines the best of both worlds? It turns out that I could.

The idea is as follows. Instead of spawning one droplet at a time and simulating it, we’ll spawn all droplets simultaneously. And not just anywhere: we spawn exactly one droplet in each grid cell. And this remains the rule throughout the simulation; each grid cell contains exactly one droplet:

This allows us to track the droplets’ positions not in a flat, 1D array, but in a 2D array aligning exactly with the height map itself. For each droplet, we track:

• the position as a floating-point number (rounding this down to the nearest integers gives the grid index)
• the speed at which it’s moving, starting out arbitrarily at 1 (notice that this isn’t really related to movement, and is more a measure of kinetic energy)
• the size (starting out at 1)
• the amount of sediment it’s carrying (starting out at 0)

Next, we compute the terrain height at each droplet’s position using bilinear interpolation, as well as the gradient of the terrain. The gradient points uphill, so we want to move the droplet in the exact opposite direction. Regardless of its speed, we always move it over a distance of exactly 1 grid step. This guarantees that it’ll either remain in its current cell, or at worst in one of the 8 neighbouring cells:

We don’t yet move the droplet’s data into the new grid cell at this point!

Next, we compute the terrain height at its new position. Subtracting the previous height gives us a height difference, i.e. how much the droplet has moved downwards (or, sometimes, upwards). We use this to calculate the droplet’s sediment carrying capacity, which is the product of the height difference, speed, size and a configurable parameter, the sediment capacity factor.

If the droplet can carry more sediment than it currently does, it picks some up from its previous location, lowering the terrain there, but never by so much that it would turn a downhill slope into an uphill one – erosion never does that. Conversely, if it carries more sediment than it can, it deposits some at its previous location, again never turning an uphill slope into a downhill one. The height change is applied to the terrain through bilinear interpolation, depending on the droplet’s position within its cell.

Next, we update the droplet’s speed depending on how far it moved downhill, and we evaporate some of its water (multiply its size by something like 0.95).

Now comes the secret sauce! We have updated each droplet’s position, but this means that some droplets (most, in fact) will have left the grid cell in which their data is stored:

This leaves some cells overcrowded, with up to 9 droplets in them! My solution is to merge these droplets into one:

We average the position and speed, weighted by the size of the original droplets, but of course we sum up the sizes and sediments without weighting.

Droplet movement also leaves some cells empty. One solution is to just spawn a new droplet with size 0 (essentially no droplet at all) in those cells. Another solution is to implement steady rainfall: in every time step, always spawn a small droplet in every cell. This is slightly more elegant because we don’t need to check for division by zero, and is the approach I ended up with.

If you want to know the details, I’ll refer you to Sebastian Lague’s implementation on GitHub. Though my data structures and outer loops are different, the equations are pretty much the same. Though I do think Sebastian’s code has some bugs…

And now the moment you’ve been waiting for: let’s see if it works!

It’s clear that the algorithm does what it should do, that it doesn’t introduce strange artifacts or numerical instabilities, and most importantly: that the result looks great! And it takes only 20 iterations, running in less than 1.5 seconds for a 2048×1024 map.

As a finishing touch, I made the erosion strength dependent on tectonic region to randomly get more or less eroded terrain depending on location. The final result:

Erosion strength is also dependent on latitude to get more hydraulic erosion towards the poles, resulting in more fjord-like structures on coastlines there (zoom in to see them).

In nature, fjords occur on terrain closer to the poles because it was once covered in ice and therefore subjected to glacial erosion. Which… I should probably not be getting into right now.