Procedure for rivers and lakes: erosion, deposition and lake formation

Generating the terrain is just half of the story: the water cycle defines the shape of the landscape we are so familiar with. Perlin noise alone lacks the features created by the flow of water in its way to the ocean. Our objective is to create something like this...

Rivers basins of the US in rainbow colours by Fejetlenfej

Rivers basins of the US in rainbow colours by Fejetlenfej


...and making the fluvial network to interact and create landscape by erosion and deposition. with a I describe here the procedure I am using based on networks: prepare your favorite graph library because we will burn it! (After some useful discussions about water tables I have started an experiment that is quite promising, some preliminary details at the end.)

1. Creating the mesh and the initial terrain

The procedure here works for any mesh (Voronoi, regular, or Delaunay as the one I will use). Delaunay triangulation is quite useful for irregular grids since it is the dual of the Voronoi one. I personally like it since I can use the Voronoi polygons for plotting any color at each point (below right, an example with 1000 points) and the Delaunay for creating the graph that defines the interaction between the grid points (below left, for the same 1000 points).


Just some basics on graphs (even if it is probably not necessary): In one hand, I have the coordinates of the mesh points that are the vertices of the Delaunay network. For dealing with the network, on the other hand, I assign an ID to every node (from 1 to N, where N is the number of points) and represent every edge by the couple of vertex IDs that are connecting. For instance, if I have a vector with all the positions there, the index of the vector works as ID and the graph is represented by another vector which elements are the couples of IDs of first neighbors. So three points (1,2,3) connected in a triangle are represented as ((1,2),(2,3),(3,1)), right?

For an easy following of the steps of the procedure, I create an artificial island composed by two volcanoes at two levels where lakes are expected in their respective calderas (for procedures creating volcanoes from tectonic activity for landscapes and planet formation check here **not finished yet!**). Note that I'm not trying to be realistic here, this is just an easy example.


2. Identifying the maximum gradient, classifying grid points, and creating the preliminary fluvial network

Gravity guides water along its path to the ocean by choosing the direction of the maximum gradient. Here, this is as easy as checking the first neighbors of every vertex and select the edge that connects it to the neighbor with the lowest elevation (strictly speaking, with the one with the biggest slope down, computed by the arctan of the difference of elevation over the horizontal distance). Direction is important here, so I save every edge of maximum gradient knowing that the direction is from the first to the second. That is, for the edge (1,2) it would mean that water flows form 1 to 2 since the elevation of 1 was higher than for 2.

You will find four cases here:

  1. Some points are in the top of a local maximum, so all its neighbors have a lower elevation. Cool, just choose that one with the biggest slope down. Two neighbors has exactly the same slope? that would be weird in random generated terrains but anyway, take both. However, for simplicity I will assume here that for every vertex the water has only one direction (ensure that you save the edges keeping the right direction). These are initial-points in the fluvial network.
  2. Some points have some neighbors with higher elevation, and others with lower. This would be the most frequent case. Just choose that one with the highest slope down as before. If it is a transit point, some of its neighbors will point to him. If no neighbor points to him because they had a better choice, that's ok no hard feelings, it will become an initial-point.
  3. Some points are in a local minimum, so there is no neighbor with lower elevation. Ok, things got interesting here, we have and end-point. Save the ID of this point but no outgoing edge comes with him. What if one neighbor has the same elevation, that is, a flat surface? Save each of them as end-points, this issue fits in the procedure below.
  4. Some points are just the ocean. Well, you can define the sea level before hand and that's ok. Classify them as ocean and that is, no need of an edge for them.

With this collection of edges we can define our preliminary fluvial network. In our example, we have one end point at the bottom of each caldera (black dots below). When you check it, something is definitely wrong: water cannot just disappear at these end-points. We need to create lakes at the end-points!


3. Merging fluvial clusters and creating lakes and seas

This is trickiest part of the procedure and it has to be done with care. But no worries, the concept is very clear. I tried many options before accepting that this was the easiest and more effective way. However, if you find an even easiest one, I will be very happy to learn! This is how I proceed:

1. I define a fluvial cluster for every end-point. As you can see in the fluvial network above, it is composed by unconnected sub-networks (the big one that includes every river that ends at the ocean, and the two that end at each end-point). Identify them using standard graph functions and libraries. I am using the next trick: I build the transpose of the adjacency matrix of the fluvial network. It is a NxN sparse matrix that I will call A. I then create a vector v with N components in such a way that all components are zero excepting the one at the position of the end-point (remember, his ID) which is 1. With the dot product A.v I got a new vector which components are zero excepting those at the first neighbors of the end-point (according to the reverse fluvial network). Thus, if I iterate with v = v + A.v I reach a limit where v doesn't change anymore and it's 1 for those nodes that belong to the cluster. You see what I'm doing here: I'm going upstream until reaching every point belonging to the river basin.

2. I include as first cluster that of the ocean. I just put together all these nodes with all the basins that end directly at the ocean. It can be done using the same trick as before, just starting in v with those points at the sea. So, at this point, we have a list of clusters, starting with the one of the ocean, with a length of number of end-points plus one. In our example, we have two end-points, that is, three clusters. The image below (1) show the three clusters: the big thing is the one of the ocean, and the other two those in each caldera.

3. I sort the edges of the original Delaunay graph according to the elevation of the highest of every couple. This will be very helpful for accelerating the next part of the procedure, I promise. Please do not underestimate this step, it's actually the key!

4. Following the clusters order, I fill them with water starting at each end-point. Now, I take the first end-point and its elevation. I fix the water level to that elevation. If I now construct the graph with only those edges below the water level, I obtain a main sub-network with all terrain below that level, and the original end-point but isolated. Indeed, he was in a minimum, so none of his first neighbors is above the water level. Next, I iteratively increase the water level by including edges one by one and following the order defined in the previous step. As I'm doing this, the graph generated by the edges below the water level will show at least two unconnected components, one of them composed by the neighbors of the original end-point. I'm showing below the example of the first volcano where you can see the two components below the water level: that of the ocean and that of the caldera. However, if I continue adding edges and increasing the water level, there will be a moment that both components will merge, as shown below right (big red arrow, you see?). This will be the path that the water will follow, and defines the final water level for this lake. This edge is actually connecting both clusters (the one of the end-point with the ocean cluster in this case), so I add it to the preliminary fluvial graph constructed before. With this I merge the clusters as show above for (2). In addition, every vertex of the cluster below the final water level will be fulled by water, so all these guys are all connected between them in the fluvial network (this will be helpful later). I repeat this operation with every cluster until all of them are merged with the ocean, as in (3) above.

Wait a second! Not every lake ends necessarily with a river connecting it with the ocean, what happen with evaporation? Yeah, right. An evaporation rate can be defined according to the climate (check my climate generation for instance). According to this, the amount of water in the caldera can be estimated to stop the rising-the-water-level procedure at the point that equilibrates with the evaporation at the surface. No prob.

There is another issue: what happen if at the bottom of a local minimum you have several end-points that actually belong to the same lake? Following the list of end-points by order, you will merge the first cluster with the second one before reaching the top water level. No prob, once you start the process with the second end-point both clusters will be merged and you will see no difference. Water level will continue growing until all points belong to the same lake at the end of the day. Be careful with the code for this part!

BTW, I promised about the flat surfaces: the end-points belonging to the flat surface will merge in the cluster creating a lake there. If they are just aligned following a line, it will form a one-dimensional lake if you want.


At the end of this part make sure that you have a list of the vertices that belong to each lake, their water level, and an updated fluvial network with the new connections between clusters and lake-vertices.

4. Obtaining the volume of flow

For the volume flow, let's use again the power of networks. With the adjacency matrix of the updated fluvial graph, I can start a path at every initial-point and go downstream counting every step I made. This is actually counting the graph-distance between the initial-point and the others following the water flow. Once it reaches the ocean, I stop. I'm assuming here that there is an equal contribution of a unit of flow for every point, so it is raining the same amount everywhere. This can be improved by the climate model defining differences in the contribution or adjusting it to the actual area of the Voronoi polygon. Fell free.

By adding the contribution from every initial-point I have the estimation of the total flow. Note that since I connected all points belonging to a lake I only need one step for entering in a lake and leaving it, meaning that as many water leaves the lake as the one that enters, in perfect equilibrium. Again, this can be modified to include higher evaporation rates, but let's assume we are in Canada or Switzerland (actually I'm in Switzerland).

I show the flow below (points at seas and lakes are represented just as dots, and I don't consider any flow at the ocean).


It is not necessary to plot every line: we can consider that there is a river only if some flow threshold is achieved.


5. Estimating erosion and deposition

Once I know the volume flow, I need to know the rule for erosion and deposition: I learned that these rates depend on the river flow and slope. The steeper, the higher the erosion; the flatter, the higher the deposition. Easy. I found useful for my example the following empirical model (but it is totally arbitrary and does not need to work properly under any circumstance, so you can be creative here):

Where Z is the increase/decrease of elevation, c is the volume flow, a is the slope from flat (a=0) to vertical (a=1) and a0 some reference value to scale it. As you see below, deposition occurs here at the coast (red) and erosion along the hillsides of the volcanoes (blue).


6. Repeat from 2 to 5

The water erosion is not a mechanism in equilibrium, so there is no convergence in this iteration. Just repeat it until you feel it is enough. I would definitely recommend to make Z small and iterate several times before making it big and use less iterations, even if it requires more computational time.

Other examples

I have below some other examples with different landscapes and forms of Z functions for erosion/deposition, just to enjoy. Changing these parameters can result in a totally different landscape, as shown there. There is a lot of room for improvements. Some important features of the procedure are, for instance: creation of fiords when erosion wins; creation of deltas when deposition wins instead; fractal-like structure of mountains and river basins after some iterations (almost independently of the initial shape); deposition along the rivers; etc. In addition, potential lakes can be filled with land instead to locate salt pans in desert areas, or landscapes with no local valleys (Perlin noise generate a lot of small valleys and a better result can be achieved if they are directly filled with land).

If you feel that some important feature are still missing with this procedure or you have any idea, just write!

Outcome for linear coastline with high erosion and low deposition.

Outcome for linear coastline with high erosion and low deposition.

Same coastline with an slightly different  Z  model. A delta is forming for the big river.

Same coastline with an slightly different Z model. A delta is forming for the big river.

Initial terrain.

Initial terrain.

Terrain after some iterations with emphasized deposition.

Terrain after some iterations with emphasized deposition.



Water Table Approach

Probably, the most realistic approach is to describe the water table itself, and consider the sponge-like behavior of land. After reading about it, I decided to work in a simple model inspired on wetting physics (I eventually have some experience on quantum liquids on porous surfaces). So here we go: I consider that water is moving between air and a porous medium under the effect of gravity. Just like exposing water to a sponge. Thus my problem is to find the surface of the water table for a given landscape. I will do it by finding the surface that minimizes the net energy function:

E = S + U - G + C

where the forces at hand are:

  1. Surface tension: It is what keeps water in the form of droplets, trying to reduce as possible the size of the exposed surface. Surface tension depends on the media, thus we have a huge tension factor when water is exposed to air, but a soft one when is in the sponge. This creates flat surfaces in lakes and seas but a more flexible one inside the sponge.
  2. Internal potential: Water molecules are bounded and they like to be together. And they are much more bounded if there is no one else in the middle, so this is actually working as a penalty for being in the sponge.
  3. Gravity. Yeah, we can not escape gravity.
  4. Capillarity. This is what makes water to climb and wet the sponge even against gravity.

There are some arbitrary parameters in the model that will change the shape of the surface. That's actually good since they can be linked with climate variables as humidity and temperature and generate different fluvial patterns depending on the climate.

I'm showing below a preliminary test in one dimension. For a given distribution of land elevation, I start with a random water table surface and iterate until converging to the one that minimizes the energy function E. Looks cool enough for trying in full 3D!

First iterations of the minimization procedure. Since the surface tension is very high when water is exposed to air, the surface at seas and lakes is formed quite fast.

First iterations of the minimization procedure. Since the surface tension is very high when water is exposed to air, the surface at seas and lakes is formed quite fast.

Subsequent iterations until convergence. Water climbs up the sponge until reach equilibrium with gravity, also affecting the final elevation of the surface of nearby lakes.

Subsequent iterations until convergence. Water climbs up the sponge until reach equilibrium with gravity, also affecting the final elevation of the surface of nearby lakes.

Same approach works in 3D indeed. To make the calculation a bit faster faster and simplifying the code for the surface tension term, I found that there is no too much difference between computing the exact value of the polygon surfaces than just using the length of the edges (anyway, if you are already using a nice library for geometry there is no such a difference).

I'm using a Hue color function for the water to clearly see the elevation of the different water bodies (red is the sea, lowest level). Only rivers ending or coming from a water body are sustainable, all the others just deposit at their end-point.


While testing values for the parameters of the energy function, I found that it can be also used for snow! cool! relax a bit the surface tension factor for water-air and voilà!