Hydraulic Erosion Algorithms in GLSL Compute

This article explores the mathematical formulation and GPU-accelerated implementation of particle-based hydraulic droplet erosion using OpenGL Shading Language (GLSL) compute shaders. It outlines the core physical equations governing droplet motion, sediment transport, deposition, and evaporation, followed by the parallel architecture and synchronization strategies required to evaluate millions of simultaneous erosion steps efficiently across a discrete digital elevation model.

Mathematical Formulation of Droplet Erosion

Particle-based hydraulic erosion models simulate discrete water droplets traveling over a continuous 2D scalar heightfield \(h(x, y)\). Each droplet maintains an internal state defined by:

Position: \(\mathbf{p} = (p_x, p_y) \in \mathbb{R}^2\)
Direction and Velocity: \(\mathbf{v} = (v_x, v_y) \in \mathbb{R}^2\)
Water Volume: \(w \in \mathbb{R}^+\)
Sediment Load: \(s \in \mathbb{R}^+\)

The simulation proceeds iteratively across discrete time steps \(\Delta t\) through the following physical phases:

1. Surface Gradient and Velocity Update

The surface gradient \(\nabla h(\mathbf{p})\) at the continuous position \(\mathbf{p}\) is obtained via bilinear interpolation of the surrounding grid cells. The terrain normal determines the acceleration vector:

\[\mathbf{a} = -\nabla h(\mathbf{p}) = -\left( \frac{\partial h}{\partial x}, \frac{\partial h}{\partial y} \right)\]

To model momentum, the velocity vector incorporates directional inertia parameter \(I \in [0, 1]\):

\[\mathbf{v}_{t+1} = \mathbf{v}_t \cdot I + \mathbf{a} \cdot (1 - I)\]

The droplet moves to its new position:

\[\mathbf{p}_{t+1} = \mathbf{p}_t + \frac{\mathbf{v}_{t+1}}{\|\mathbf{v}_{t+1}\|} \cdot \Delta t\]

2. Sediment Capacity Calculation

The maximum sediment capacity \(C\) a droplet can carry depends on its speed, water volume, and the local slope \(\Delta h = h(\mathbf{p}_{t+1}) - h(\mathbf{p}_t)\):

\[C = \max(-\Delta h, S_{\min}) \times \|\mathbf{v}\| \times w \times K_c\]

where \(S_{\min}\) is a minimum slope threshold to prevent sediment stagnation on flat regions, and \(K_c\) is a sediment capacity scaling constant.

3. Erosion and Deposition

The sediment differential determines whether erosion or deposition occurs:

Deposition (\(s > C\)): If the droplet carries more sediment than \(C\), a fraction proportional to the deposition rate \(K_d\) is placed back onto the terrain:

\[\Delta s = (s - C) \times K_d\] \[s_{t+1} = s_t - \Delta s\] \[h(\mathbf{p}) \leftarrow h(\mathbf{p}) + \Delta s\]

Erosion (\(s < C\)): If the droplet carries less sediment than \(C\) and is traveling downhill (\(\Delta h < 0\)), terrain height is removed up to the capacity deficit, regulated by the erosion rate \(K_e\):

\[\Delta s = \min((C - s) \times K_e, -\Delta h)\] \[s_{t+1} = s_t + \Delta s\] \[h(\mathbf{p}) \leftarrow h(\mathbf{p}) - \Delta s\]

The mass change \(\Delta s\) is distributed among or subtracted from the four neighboring heightfield grid cells using bilinear weighting factors.

4. Evaporation and Entrainment

At the end of each iteration, the droplet loses water according to an evaporation rate \(K_e \in (0, 1)\), and its speed decreases due to friction:

\[w_{t+1} = w_t \times (1 - K_e)\] \[\mathbf{v}_{t+1} = \mathbf{v}_{t+1} \times (1 - \text{friction})\]

GLSL Compute Shader Implementation

Implementing this algorithm in GLSL compute shaders introduces two execution strategies: droplet-centric compute and grid-centric (shallow water) simulation. The droplet-centric approach maps each compute thread to an individual droplet life cycle.

Managing Concurrent Memory Access

Because multiple droplet threads can alter the same heightfield cells simultaneously, race conditions occur when modifying the elevation buffer. In GLSL, this is handled using an integer-based fixed-point representation stored in a Shader Storage Buffer Object (SSBO) with atomic operations:

\[\text{FixedPointInt} = \text{int}(h \times 1000.0)\]

This enables deterministic read-modify-write cycles using atomicAdd:

void modifyHeight(ivec2 coord, float deltaH) {
    int fixedDelta = int(deltaH * 1000.0);
    int index = coord.y * MAP_WIDTH + coord.x;
    atomicAdd(heightMapBuffer[index], fixedDelta);
}

Core Droplet Execution Loop

The compute shader manages the path and terrain adjustments within a single kernel dispatch:

#version 450 core

layout(local_size_x = 64, local_size_y = 1, local_size_z = 1) in;

layout(std430, binding = 0) buffer Heightfield {
    int heightData[]; // Scaled fixed-point integers
};

uniform int MAP_WIDTH;
uniform int MAP_HEIGHT;
uniform int MAX_LIFETIME;
uniform float INERTIA;
uniform float CAPACITY_FACTOR;
uniform float EROSION_RATE;
uniform float DEPOSITION_RATE;
uniform float EVAPORATION_RATE;

// Pseudo-random position generation per invocation
vec2 getRandomPosition(uint id);

// Bilinear gradient and height extraction
vec3 getSurfaceNormalAndHeight(vec2 pos);

void main() {
    uint id = gl_GlobalInvocationID.x;
    vec2 pos = getRandomPosition(id);
    vec2 dir = vec2(0.0);
    
    float water = 1.0;
    float sediment = 0.0;
    float speed = 1.0;

    for (int step = 0; step < MAX_LIFETIME; ++step) {
        ivec2 iPos = ivec2(pos);
        if (iPos.x < 1 || iPos.x >= MAP_WIDTH - 2 || 
            iPos.y < 1 || iPos.y >= MAP_HEIGHT - 2) {
            break;
        }

        vec3 surfaceData = getSurfaceNormalAndHeight(pos);
        vec2 grad = surfaceData.xy;
        float hOld = surfaceData.z;

        dir = (dir * INERTIA) - (grad * (1.0 - INERTIA));
        float len = length(dir);
        if (len == 0.0) break;
        dir /= len;

        vec2 newPos = pos + dir;
        float hNew = getSurfaceNormalAndHeight(newPos).z;
        float deltaH = hNew - hOld;

        float capacity = max(-deltaH, 0.01) * speed * water * CAPACITY_FACTOR;

        if (sediment > capacity || deltaH > 0.0) {
            float depositAmount = (deltaH > 0.0) ? min(deltaH, sediment) : (sediment - capacity) * DEPOSITION_RATE;
            sediment -= depositAmount;

            // Distribute deposit over four neighbor cells via bilinear weights
            vec2 f = fract(pos);
            modifyHeight(iPos, depositAmount * (1.0 - f.x) * (1.0 - f.y));
            modifyHeight(iPos + ivec2(1, 0), depositAmount * f.x * (1.0 - f.y));
            modifyHeight(iPos + ivec2(0, 1), depositAmount * (1.0 - f.x) * f.y);
            modifyHeight(iPos + ivec2(1, 1), depositAmount * f.x * f.y);
        } else {
            float erodeAmount = min((capacity - sediment) * EROSION_RATE, -deltaH);
            sediment += erodeAmount;

            vec2 f = fract(pos);
            modifyHeight(iPos, -erodeAmount * (1.0 - f.x) * (1.0 - f.y));
            modifyHeight(iPos + ivec2(1, 0), -erodeAmount * f.x * (1.0 - f.y));
            modifyHeight(iPos + ivec2(0, 1), -erodeAmount * (1.0 - f.x) * f.y);
            modifyHeight(iPos + ivec2(1, 1), -erodeAmount * f.x * f.y);
        }

        speed = sqrt(max(0.0, speed * speed + deltaH * 9.81));
        water *= (1.0 - EVAPORATION_RATE);
        pos = newPos;

        if (water < 0.01) break;
    }
}

Performance Considerations

Memory Coalescing: Random droplet movement creates high divergence in memory access. Thread workgroups run more efficiently when spawned within spatially localized bounding boxes rather than across random disjoint intervals.
Atomic Contention: If hundreds of droplets converge into the same drainage channels, atomic serialization limits throughput. Splitting iteration steps into staggered dispatch passes reduces concurrent memory collisions on critical path height nodes.