Simulating Terrain Erosion with GLSL Compute Shaders

Simulating sediment erosion and deposition on terrain heightmaps using GLSL compute shaders allows developers to leverage the massive parallelism of modern GPUs to generate realistic, geologically plausible landscapes in real time. By representing the terrain as a two-dimensional grid of height and sediment values, compute shaders can iteratively model the physical interactions of rainfall, surface water flow, soil dissolution, sediment transport, and deposition across thousands of virtual landscape cells simultaneously.

Heightmap Representation and Texture Buffers

In a GLSL-based erosion framework, terrain is typically stored in floating-point 2D textures or Shader Storage Buffer Objects (SSBOs). A standard grid-based Eulerian approach utilizes multiple data channels across one or more textures:

  • Terrain Height (\(H\)): The base bedrock and consolidated rock elevation.
  • Water Level (\(W\)): The surface depth of the water flowing over the terrain.
  • Suspended Sediment (\(S\)): The amount of eroded soil currently carried by water.
  • Outflow Flux (\(F\)): A four-component vector representing water flowing to the four orthogonal neighbor cells (North, South, East, West).
  • Velocity Vector (\(V\)): The directional movement and speed of surface water.

Compute shaders read and write to these layers using layout(rgba32f) uniform image2D qualifiers alongside imageLoad() and imageStore() operations.

The Multi-Pass Compute Pipeline

Because fluid and soil mechanics depend on neighboring texels, the simulation cannot execute safely in a single pass without introducing race conditions. Instead, an erosion cycle is broken down into sequential compute shader dispatches, separated by host-side memory barriers (glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT)):

1. Water Influx

A compute pass simulates precipitation by incrementally adding a specified rainfall amount to the water level texture:

layout(local_size_x = 16, local_size_y = 16) in;

layout(binding = 0, r32f) uniform image2D waterTexture;
uniform float u_RainRate;
uniform float u_DeltaTime;

void main() {
    ivec2 coord = ivec2(gl_GlobalInvocationID.xy);
    float water = imageLoad(waterTexture, coord).r;
    water += u_RainRate * u_DeltaTime;
    imageStore(waterTexture, coord, vec4(water, 0.0, 0.0, 0.0));
}

2. Flow and Flux Calculation

Water movement is typically calculated using the pipe model. The shader evaluates hydrostatic pressure differences between the current cell and its four neighbors: \[\Delta h = (H_{current} + W_{current}) - (H_{neighbor} + W_{neighbor})\] If \(\Delta h > 0\), water flows outward. The new flux to each neighbor is calculated, scaled by gravity and cross-sectional area. If the sum of all outward fluxes exceeds the total water volume available in the cell, the output fluxes are scaled down uniformly to prevent negative water depths.

3. Erosion and Deposition

Velocity is derived from the net flow of water across the cell boundaries. The compute shader calculates the local slope and fluid speed, which together define the sediment transport capacity (\(C\)): \[C = K_c \cdot \sin(\theta) \cdot \|V\|\] Where \(K_c\) is the sediment capacity constant, \(\theta\) is the surface slope angle, and \(\|V\|\) is the flow velocity.

  • If Suspended Sediment (\(S\)) < Capacity (\(C\)): The fluid is under-saturated. The shader erodes bedrock by moving a fraction of material from the heightmap layer into the sediment layer: \[Erosion = K_e \cdot (C - S)\]
  • If Suspended Sediment (\(S\)) > Capacity (\(C\)): The fluid is over-saturated. Excess sediment drops out of the water and permanently hardens back onto the heightmap: \[Deposition = K_d \cdot (S - C)\]

4. Sediment and Water Advection

Sediment must physically move along with the surface velocity. Instead of pushing sediment forward (which causes race conditions), the shader performs semi-Lagrangian back-tracing:

vec2 previousPos = vec2(coord) - velocity * u_DeltaTime;
float transportedSediment = texture(sedimentSampler, previousPos / textureSize).r;

Sampling the sediment texture using hardware-accelerated bilinear interpolation produces smooth advection artifacts rather than hard, blocky grid edges.

5. Evaporation

Finally, a fraction of the water volume is evaporated based on an evaporation constant (\(K_{evap}\)): \[W = W \cdot (1.0 - K_{evap} \cdot \Delta t)\] As water depth approaches zero, any remaining sediment carried by that water is fully deposited back onto the terrain surface, forming alluvial fans at the bases of valleys and river outlets.

Synchronization and Memory Management

Compute shaders execute work in localized workgroups (typically 16x16 or 32x32 threads). Because neighbor reads are critical for gradient and flux equations, developers commonly utilize ping-pong texturing—alternating between two textures for reading and writing on successive dispatches—to eliminate write hazards. For systems leveraging OpenGL 4.3+ or Vulkan, this compute pipeline processes millions of grid cells across hundreds of simulated years in a matter of seconds.