☂️ 💏 💂🏼 Particle based surface erosion simulation 💖 🦉 🌌

Note: the full source code of the project, as well as explanations about its use and reading, can be found on Github [ here ].

I took a break from my master's thesis to work on what I had been putting off for a long time: improved terrain generation for my Territory project . A simple way to implement it is hydraulic erosion, which is why I created it!

For a one-day software puzzle, it worked pretty well, and was not as complicated as I expected. Results are quickly generated, have a physical meaning, and look amazing.

In this article, I will talk about my simple C ++ implementation of a particle-based square mesh hydraulic erosion system. I will explain all the physical justifications underlying the implementation, and talk about mathematics. The code is extremely simple (only about 20 lines for the mathematics of erosion) and quick to implement, so I recommend it to anyone who wants to increase the realism of their terrain.

The results are rendered using a stripped down version of my Homebrew OpenGl Engine , which I modified to render a 2D array of points as a height map. A stripped-down version of the engine is much easier to understand if you are interested in learning OpenGL in C ++.

Implementation

Inspired by many sources of hydraulic erosion, I decided that it would be most logical to use particles.

Particle-based erosion is very simple:

We create a particle at a random point on the surface.
It moves / slides along the surface using standard classical mechanics (we will talk about them below)
We carry out the transfer of matter / sediment between the surface and the particle (this is also explained below)
We evaporate part of the particle
If the particle is outside the map or too small, then destroy it
Repeat the process with the desired number of particles.

Note: system parameters are explained in the corresponding section. The most important is the time step factor dt, which proportionally scales all parameters. This allows us to increase the frequency of the simulation (at the cost of increasing noise) without changing the relative scale of the parameters. This can be seen in the code below.

Particles

To do this, I created a simple particle structure containing all the properties I need:

struct Particle{
  //Construct particle at position _pos
  Particle(glm::vec2 _pos){ pos = _pos; }

  glm::vec2 pos;
  glm::vec2 speed = glm::vec2(0.0); //Initialize to 0

  float volume = 1.0;   //Total particle volume
  float sediment = 0.0; //Fraction of volume that is sediment!
};

Note: in case you are not familiar with the GLM library: I use it to perform vector operations.

A particle has a position and speed that determines how it moves. In addition, it has a volume and a fraction that determines how much of the volume is sedimentary rocks.

Movement: classic mechanics

Particle movement on the surface is simulated using classical mechanics . In short, the position x of the particle is changed by the velocity v , changed by the acceleration a .

Note: bold letters indicate that the value is a vector.

We also know that force is equal to mass times acceleration:

The particle experiences downward acceleration caused by gravity, but it is on the surface, which makes downward acceleration impossible. Therefore, instead of this, the particle is subjected to a force F directed along the surface and proportional to the normal to the surface.

Therefore, we can say that the acceleration a is proportional to the normal vector of the surface n divided by the mass of the particle.

where k is the proportionality constant, and m is the particle mass. If the mass is equal to the volume multiplied by the density, then we obtain the complete system of particle motion using classical mechanics:

//... particle "drop" was spawned above at random position

glm::ivec2 ipos = drop.pos; //Floored Droplet Initial Position
glm::vec3 n = surfaceNormal(ipos.x, ipos.y);  //Surface Normal

//Accelerate particle using classical mechanics
drop.speed += dt*glm::vec2(n.x, n.z)/(drop.volume*density);
drop.pos   += dt*drop.speed;
drop.speed *= (1.0-dt*friction);  //Friction Factor

//...

Note: speed after particle motion is reduced by the friction vector. Notice that the time step factor is included here. Particles have their own inertia, proportional to their density due to the fact that we simulate their motion using acceleration.

The process of formation of sedimentary rocks: mass transfer

The process of the formation of sedimentary rocks physically occurs as the transfer of sedimentary rocks from the earth to the particle and back to the particle location (" mass transfer ").

In chemical technology, mass transfer (i.e. the change in mass / sediment over time) between two phases (in this case, the surface of the earth and a drop) is usually described using mass transfer coefficients .

Mass transfer is proportional to the difference between concentration c and equilibrium concentration c_eq:

where k is the constant of proportionality (mass transfer coefficient). This difference between equilibrium and actual concentration is often called the “driving force."

If the equilibrium concentration is higher than the current, then the particle absorbs sedimentary rocks. If lower, then loses them. If they are equal, then no changes will occur.

The mass transfer coefficient can be interpreted in various ways:

As the frequency of the transition between phases (here it is “deposition rate”)
The rate at which droplet concentration tends to equilibrium concentration

The system is completely dependent on determining the equilibrium concentration. Depending on the definition, the system will demonstrate different dynamics of sedimentary rock deposition. In my implementation, the equilibrium concentration is higher if we move down, and if we move faster, and also it is proportional to the volume of the particle:

//...

//Compute Equilibrium Sediment Content
float c_eq = drop.volume*glm::length(drop.speed)*(heightmap[ipos.x][ipos.y]-heightmap[(int)drop.pos.x][(int)drop.pos.y]);

if(c_eq < 0.0) c_eq = 0.0;

//Compute Capacity Difference ("Driving Force")
float cdiff = c_eq - drop.sediment;

//Perform the Mass Transfer!
drop.sediment += dt*depositionRate*cdiff;
heightmap[ipos.x][ipos.y] -= dt*drop.volume*depositionRate*cdiff;

//...

Note: the change in concentration inside the particle is completely described by the mass transfer equation. The change in the height map is additionally multiplied by the volume of particles, because we change it proportionally not to concentration, but to mass (concentration is multiplied by volume).

Other aspects

The height map is initialized by a multi-layer Perlin noise with a random seed.

At the end of each time step, the particle loses a little mass in accordance with the evaporation rate:

//...

drop.volume *= (1.0-dt*evapRate);

//...

This process is repeated for thousands of particles created in random places and simulated separately (in my case, the calculations are performed sequentially in the CPU).

My default implementation includes a good set of parameters. After 200 thousand particles, erosion looks very good.

results

The finished code for the erosion process is approximately 20 lines without comment.

Here is a “before and after” comparison for 10 samples. The simulation generates very beautiful ridges on elevations, sedimentary rocks are deposited on the sides of some ridges, which leads to the creation of beautiful plateaus.

A selection of ten before-and-after comparisons. The shader I wrote receives two colors at the input. I also implemented shadow mapping, distance-dependent fog and Phong shading.

Here are ten more (different) samples just with the results:

Ten more sample results. They differ from the previous ones, even if some formations look similar.

By choosing seed differently during map initialization, you can create controlled different results.

Simulation time

Simulation time is directly related to the lifetime of the particles and the number of simulated particles.

Several factors affect the life of particles, and it can vary greatly for each particle, because they are created in random places:

Friction and inertia: particle velocity
Grid Size: The probability that particles will fall out of the map.
Evaporation rate: particle extinction rate

With the default parameters, a simulation of an individual particle requires 10-100 milliseconds, which results in 10-20 seconds for a simulation of 200,000 particles.

You can increase the degree of erosion and reduce the life of the particles by increasing the time step without changing the number of simulated particles. This can make the simulation more “noisy”, and if the time step is too large, there is a chance that the simulation will fail.

The simulation itself in the engine incurs additional costs for re-creating the surface mesh for rendering.

You can optimize the code by reducing the cost of re-creating the mesh or by increasing the speed of one time step for an individual particle.

Work for the future

Soon I will insert this code into my Territory project as the basis for generating a height map in the elevation generator.

I already wrote a framework for simplified simulation of fluid dynamics in climate systems, but it is not yet ready for publication. This system gets climate patterns from a height map. This will be the topic for a future post on physically accurate procedural climate systems!

Using a simulated climate system, it will be possible in the future to sample particles from the distribution (for example, in places where it rains) instead of distributing them evenly on the map. They will be ideally combined with the geology and climate of the relief.

In addition, after adding different types of geological formations, stones can have a different degree of solubility. This will directly affect the mass transfer coefficient (deposition rate), creating a different erosion rate.

This system is unable to simulate a real fluid and river flow, but it can potentially be adapted to such tasks. I’ll think it over and maybe in the future I will release a continuation of the post.

Particle based surface erosion simulation