Fluid Rendering in Games — Heightmap Method

1. Wave Equation for the Heightmap

The fluid surface — a thin layer with small perturbations of amplitude h(x,y,t), where h is the deviation from equilibrium. In the linear approximation h satisfies the wave equation:

This is the linear wave equation — an incompressible surface wave in the long-wave limit (λ ≫ d). It produces a ripple effect: circles spread at constant speed, reflect from walls, interfere, and decay through damping.

What It Does NOT Describe

Nonlinear effects (breaking waves, spindrift, bore waves), vertical particle motion (important for large amplitudes), full Navier-Stokes dynamics. But for interactive games — sufficient.

2. Explicit Finite-Difference Scheme

We discretise in space (steps Δx, Δy) and time (step Δt). Notation: h^n_{i,j} = h(i·Δx, j·Δy, n·Δt). Explicit "leap-frog" scheme:

Violating the Courant-Friedrichs-Lewy condition causes the scheme to diverge — numerical oscillations grow. In games typically: Δt = 1/60 s, Δx = tile size. Maximum wave speed c_max ≤ Δx·60/√2.

Damping

Without damping waves reflect forever — unrealistic. Multiply by a damping coefficient d < 1 after each step:

3. Cellular Automaton Approximation

There is an even simpler but effective version — a cellular automaton approach popularised in Lode Vandevenne's blog. We store only two textures: current and previous height. Update rule:

// Fragment shader — one wave cellular automaton step
uniform sampler2D uCurr;   // current state h^n
uniform sampler2D uPrev;   // previous state h^{n-1}
uniform vec2 uTexelSize;   // 1/N for each direction
uniform float uSpeed;     // (c*dt/dx)^2
uniform float uDamping;   // 0.99
varying vec2 vUV;

void main() {
  float c = texture2D(uCurr, vUV).r;
  float p = texture2D(uPrev, vUV).r;
  // Neighbours
  float n = texture2D(uCurr, vUV + vec2(0, uTexelSize.y)).r;
  float s = texture2D(uCurr, vUV - vec2(0, uTexelSize.y)).r;
  float e = texture2D(uCurr, vUV + vec2(uTexelSize.x, 0)).r;
  float w = texture2D(uCurr, vUV - vec2(uTexelSize.x, 0)).r;

  float laplacian = n + s + e + w - 4.0 * c;
  float next = (2.0 * c - p + uSpeed * laplacian) * uDamping;
  gl_FragColor = vec4(next, 0., 0., 1.);
}

This shader runs full-screen, ping-pong between two FBOs — identical to a reaction-diffusion simulation. Operations count: O(N²) per frame on GPU — even a 512×512 grid runs comfortably at 60 FPS.

Interaction: Throwing a Stone

function splashAt(gl, fbo, x, y, radius, amplitude) {
  // Draw a Gaussian "bump" at the contact point
  // Then redraw into the current FBO with additive blending
  gl.bindFramebuffer(gl.FRAMEBUFFER, fbo);
  gl.enable(gl.BLEND);
  gl.blendFunc(gl.ONE, gl.ONE);
  gl.useProgram(splatProg);
  gl.uniform2f(loc('uCenter'), x, y);
  gl.uniform1f(loc('uRadius'), radius);
  gl.uniform1f(loc('uAmp'), amplitude);
  drawFullscreenQuad(gl);
  gl.disable(gl.BLEND);
}

4. Foam Map via Velocity Accumulation

Foam forms where water moves fast or where wave crests break. The heightmap model has no real tumbling, but we can approximate it:

// In the foam update shader:
float vel = abs(next - p);               // local height change
float newFoam = clamp((vel - 0.005) * 20.0, 0.0, 1.0);
float oldFoam = texture2D(uFoam, vUV).r;
float foam = mix(oldFoam, newFoam, 0.1) * 0.96; // decay
gl_FragColor = vec4(next, foam, 0.0, 1.0); // r=height, g=foam

We pack height into the R channel and foam into the G channel — one texture instead of two. In the final render shader, foam is used as a mask for white foam on top of the water colour.

5. Normals from the Heightmap

The normal to the surface z = h(x,y) at point (x,y) is determined by the gradient: n = normalize(−∂h/∂x, −∂h/∂y, 1).

// In vertex or fragment shader:
float ts = uTexelSize.x;
float hR = texture2D(uHeight, uv + vec2( ts, 0.)).r;
float hL = texture2D(uHeight, uv + vec2(-ts, 0.)).r;
float hU = texture2D(uHeight, uv + vec2(0.,  ts)).r;
float hD = texture2D(uHeight, uv + vec2(0., -ts)).r;
vec3 normal = normalize(vec3(hL - hR, hD - hU, 2.0 * ts * uHeightScale));

The uHeightScale value matters: too large gives a "metallic" look with harsh highlights, too small — flat water. Typically 0.5–2.0 depending on grid size and wave amplitudes.

6. Fresnel: Mirror vs Transparency

Real water at a shallow viewing angle is nearly a mirror; at a large angle — transparent. This is the Fresnel effect (familiar from the caustics article):

F = 1 → reflect environment map (sky, clouds). F = 0 → look through the water (bottom, underwater objects).

float cosTheta = max(0.0, dot(normal, viewDir));
float F0 = 0.02;
float fresnel = F0 + (1.0 - F0) * pow(1.0 - cosTheta, 5.0);
// Final colour = blend of reflected and refracted
vec3 color = mix(refractionColor, reflectionColor, fresnel);

7. Beer's Law and Refraction UV Offset

Beer's Law (colour absorption)

Colour passing through water is absorbed differently for R, G, B components. Beer-Lambert law:

float depth = max(0.0, uWaterDepth - refractionZ);
vec3 absorption = exp(-vec3(0.45, 0.07, 0.03) * depth);
vec3 refractionColor = bottomColor * absorption;
// + deep water tint (bottom invisible at large depths)
refractionColor = mix(deepWaterColor, refractionColor, exp(-depth * 0.5));

Refraction UV Offset (screen-space)

Instead of true refracted ray tracing: offset the UV look-up into the screen-space bottom texture by a vector proportional to the projection of the normal onto the XY plane:

float refractStrength = 0.02;  // refraction strength
vec2 refractOffset = normal.xy * refractStrength;
// Offset UV for sampling the "underwater" texture
vec2 refractUV = screenUV + refractOffset;
vec3 bottomColor = texture2D(uBottom, refractUV).rgb;

This is not physically exact, but looks convincing and is cheap — one additional texture lookup per pixel.

8. Full GLSL Water Rendering Shader

// Water surface fragment shader
precision highp float;
uniform sampler2D uHeightFoam;  // R=height, G=foam
uniform sampler2D uEnvMap;      // cubemap or spherical sky texture
uniform sampler2D uBottom;      // bottom texture
uniform vec3  uSunDir;          // direction to sun
uniform vec3  uViewPos;
uniform float uWaterDepth;
uniform float uHeightScale;
uniform vec2  uTexelSize;
varying vec2  vUV;
varying vec3  vWorldPos;

void main() {
  float ts = uTexelSize.x;
  float hR = texture2D(uHeightFoam, vUV+vec2( ts,0)).r;
  float hL = texture2D(uHeightFoam, vUV-vec2( ts,0)).r;
  float hU = texture2D(uHeightFoam, vUV+vec2(0, ts)).r;
  float hD = texture2D(uHeightFoam, vUV-vec2(0, ts)).r;
  float foam = texture2D(uHeightFoam, vUV).g;
  vec3 normal = normalize(vec3(hL-hR, hD-hU, 2.0*ts*uHeightScale));

  vec3 viewDir = normalize(uViewPos - vWorldPos);
  float cosV = max(0.0, dot(normal, viewDir));

  // Fresnel
  float fresnel = 0.02 + 0.98 * pow(1.0 - cosV, 5.0);

  // Reflection (environment map)
  vec3 reflDir = reflect(-viewDir, normal);
  vec3 reflection = texture2D(uEnvMap, reflDir.xy * 0.5 + 0.5).rgb;

  // Refraction (screen-space approximation)
  vec2 refractUV = vUV + normal.xy * 0.025;
  vec3 bottomCol = texture2D(uBottom, clamp(refractUV, 0.001, 0.999)).rgb;

  // Beer's Law absorption
  float depth = uWaterDepth;
  vec3 absorbed = bottomCol * exp(-vec3(0.45, 0.07, 0.03) * depth);
  vec3 deepTint  = vec3(0.02, 0.08, 0.15); // deep ocean
  vec3 refraction = mix(deepTint, absorbed, exp(-depth * 0.3));

  // Specular (Blinn-Phong for sun highlight)
  vec3 halfDir = normalize(viewDir + uSunDir);
  float spec = pow(max(0.0, dot(normal, halfDir)), 256.0) * 2.0;

  // Final combination
  vec3 color = mix(refraction, reflection, fresnel) + spec;
  // Foam — always on top
  color = mix(color, vec3(1.0), foam * 0.8);

  gl_FragColor = vec4(pow(color, vec3(1.0/2.2)), 1.0); // gamma
}

Comparison of Approaches for Game Fluid

🌊 Open Ocean Simulation →