Procedural Terrain Generation with Erosion Simulation

Raw Perlin noise makes bumpy, uniform hills — not the ridgelines, valleys, and river deltas of real landscapes. This tutorial builds a complete terrain pipeline: layered Simplex noise for the base heightfield, domain warping for organic distortion, a droplet-based hydraulic erosion simulation for realistic valleys and sediment fans, thermal erosion for scree slopes, and a normal-mapped mesh ready to render in Three.js.

1. The Heightfield Representation

A heightfield is a 2D grid of floats height[x][z] representing elevation. For a resolution of N×N, a Float32Array of length N² is dramatically cheaper to simulate than a full voxel volume, and it maps directly onto a PlaneGeometry in Three.js by displacing each vertex along Y.

class Heightfield {
  constructor(size) {
    this.size = size; // grid is size × size
    this.data = new Float32Array(size * size);
  }
  get(x, z) {
    x = Math.min(Math.max(x, 0), this.size - 1);
    z = Math.min(Math.max(z, 0), this.size - 1);
    return this.data[z * this.size + x];
  }
  set(x, z, h) {
    this.data[z * this.size + x] = h;
  }
  // Bilinear sample at fractional coordinates — used heavily by erosion
  sample(x, z) {
    const x0 = Math.floor(x), z0 = Math.floor(z);
    const fx = x - x0, fz = z - z0;
    const h00 = this.get(x0, z0),     h10 = this.get(x0+1, z0);
    const h01 = this.get(x0, z0+1), h11 = this.get(x0+1, z0+1);
    const a = h00 * (1-fx) + h10 * fx;
    const b = h01 * (1-fx) + h11 * fx;
    return a * (1-fz) + b * fz;
  }
}

2. Layered Noise: Fractal Brownian Motion

A single octave of Simplex noise produces smooth, rounded hills with no fine detail. Fractal Brownian Motion (fBm) sums multiple octaves of noise at increasing frequency and decreasing amplitude, adding large landmasses, medium hills, and small-scale roughness in one pass:

fBm(x, z) = Σ (k=0..octaves-1) amplitude_k \cdot noise(x \cdot freq_k, z \cdot freq_k) freq_k = lacunarity^k (typically 2.0 — each octave doubles frequency) amplitude_k = persistence^k (typically 0.5 — each octave halves contribution)

function fbm(x, z, {
  octaves = 6, lacunarity = 2.0, persistence = 0.5, scale = 0.01
} = {}) {
  let total = 0, amplitude = 1, frequency = 1, maxAmp = 0;
  for (let o = 0; o < octaves; o++) {
    total += simplex2(x * scale * frequency, z * scale * frequency) * amplitude;
    maxAmp += amplitude;
    amplitude *= persistence;
    frequency *= lacunarity;
  }
  return total / maxAmp; // normalized to roughly [-1, 1]
}

Ridged noise variant: for mountain ranges, take 1 - Math.abs(noise) per octave instead of raw noise. This folds valleys into sharp ridgelines and is the standard trick behind most "mountain generator" demos.

3. Domain Warping

Plain fBm still looks grid-aligned and repetitive at a glance. Domain warping distorts the input coordinates with another noise field before sampling, breaking up straight ridgelines into organic, wandering shapes — the technique popularized by Inigo Quilez:

function warpedHeight(x, z) {
  // Sample a low-frequency noise field to get a distortion vector
  const warpX = fbm(x + 5.2, z + 1.3, { octaves: 3, scale: 0.02 });
  const warpZ = fbm(x + 9.1, z + 7.4, { octaves: 3, scale: 0.02 });
  const strength = 40; // world-space warp amplitude

  return fbm(
    x + warpX * strength,
    z + warpZ * strength,
    { octaves: 6, scale: 0.01 }
  );
}

4. Hydraulic Erosion: Droplet Simulation

Real terrain is carved by millions of years of rainfall. The fastest practical approximation is particle-based hydraulic erosion (Hans Beyer's droplet algorithm, also used in World Machine and Unity terrain tools): simulate a raindrop that flows downhill, picks up sediment where it moves fast, and deposits sediment where it slows down.

Per droplet, per step: 1. Compute gradient \nablah at droplet position (bilinear) 2. dir = dir\cdotinertia - \nablah\cdot(1-inertia), normalize 3. Move droplet by dir \cdot stepSize 4. Δh = height(new) - height(old) 5. capacity = max(-Δh, minSlope) \cdot speed \cdot water \cdot sedimentCapacity 6. if sediment > capacity \to deposit (sediment - capacity) \cdot deposition if sediment \leq capacity \to erode min((capacity - sediment) \cdot erosion, -Δh) 7. speed = sqrt(speed² + Δh\cdotgravity) [if downhill] 8. water *= (1 - evaporation)

function simulateDroplet(field, x, z, params) {
  let pos = { x, z };
  let dir = { x: 0, z: 0 };
  let speed = 1, water = 1, sediment = 0;

  for (let step = 0; step < params.maxSteps; step++) {
    const grad = gradient(field, pos.x, pos.z);
    dir.x = dir.x * params.inertia - grad.x * (1 - params.inertia);
    dir.z = dir.z * params.inertia - grad.z * (1 - params.inertia);
    const len = Math.hypot(dir.x, dir.z) || 1;
    dir.x /= len; dir.z /= len;

    const oldHeight = field.sample(pos.x, pos.z);
    pos.x += dir.x; pos.z += dir.z;
    if (!inBounds(field, pos)) break;

    const newHeight = field.sample(pos.x, pos.z);
    const deltaH = newHeight - oldHeight;

    const capacity = Math.max(-deltaH, params.minSlope)
      * speed * water * params.sedimentCapacity;

    if (deltaH > 0 || sediment > capacity) {
      // uphill or overloaded → deposit
      const deposit = deltaH > 0
        ? Math.min(deltaH, sediment)
        : (sediment - capacity) * params.depositionRate;
      sediment -= deposit;
      addHeight(field, pos.x, pos.z, deposit);
    } else {
      // downhill and under capacity → erode
      const erode = Math.min((capacity - sediment) * params.erosionRate, -deltaH);
      sediment += erode;
      addHeight(field, pos.x, pos.z, -erode);
    }

    speed = Math.sqrt(Math.max(0, speed * speed + deltaH * params.gravity));
    water *= (1 - params.evaporation);
    if (water < 0.01) break;
  }
}

// Run 50,000–200,000 droplets from random surface points for a 512² map
function erode(field, dropletCount, params) {
  for (let i = 0; i < dropletCount; i++) {
    const x = Math.random() * field.size;
    const z = Math.random() * field.size;
    simulateDroplet(field, x, z, params);
  }
}

Deposition uses bilinear splatting: addHeight distributes the deposited/eroded amount across the 4 surrounding grid cells weighted by the same bilinear weights used to sample — otherwise the droplet snaps to integer cells and leaves visible grid artifacts.

5. Thermal Erosion: Talus Collapse

Hydraulic erosion alone leaves slopes steeper than any real material can hold. Thermal erosion models gravity moving loose material downhill whenever a slope exceeds the angle of repose (talus angle) — around 30-40° for dry soil and rock scree:

function thermalErosion(field, talusAngle = 0.6, amount = 0.5) {
  const next = field.data.slice(); // double-buffer: read old, write new

  for (let z = 1; z < field.size - 1; z++) {
    for (let x = 1; x < field.size - 1; x++) {
      const h = field.get(x, z);
      let maxDrop = 0, targetX = x, targetZ = z;

      // Find the steepest downhill neighbor among 8 neighbors
      for (const [dx, dz] of NEIGHBORS_8) {
        const nh = field.get(x + dx, z + dz);
        const drop = h - nh;
        if (drop > maxDrop) { maxDrop = drop; targetX = x + dx; targetZ = z + dz; }
      }

      if (maxDrop > talusAngle) {
        // Move half the excess above the talus angle to the lowest neighbor
        const transfer = (maxDrop - talusAngle) * amount * 0.5;
        next[z * field.size + x] -= transfer;
        next[targetZ * field.size + targetX] += transfer;
      }
    }
  }
  field.data = next;
}

Interleave the passes: run a handful of thermal erosion iterations after every batch of ~1,000 droplets. Alternating hydraulic and thermal erosion produces the characteristic mix of smooth valleys and sharp, stable ridgelines seen in real mountain ranges.

6. Flow Accumulation and Rivers

To render visible rivers, compute a flow accumulation map: process cells from highest to lowest elevation, and at each cell push its water quantity to its single steepest downhill neighbor (the D8 algorithm). Cells with high accumulated flow are rivers; the rest are hillside.

function flowAccumulation(field) {
  const n = field.size * field.size;
  const flow = new Float32Array(n).fill(1); // every cell starts with 1 unit of rain
  const order = [...Array(n).keys()].sort(
    (a, b) => field.data[b] - field.data[a] // highest first
  );

  for (const i of order) {
    const x = i % field.size, z = Math.floor(i / field.size);
    const h = field.data[i];
    let steepest = -Infinity, targetIdx = -1;

    for (const [dx, dz] of NEIGHBORS_8) {
      const nx = x + dx, nz = z + dz;
      if (!inBoundsXZ(field, nx, nz)) continue;
      const drop = (h - field.get(nx, nz)) / Math.hypot(dx, dz);
      if (drop > steepest) { steepest = drop; targetIdx = nz * field.size + nx; }
    }
    if (targetIdx >= 0 && steepest > 0) flow[targetIdx] += flow[i];
  }
  return flow; // log(flow) > threshold ⇒ paint as river in the terrain texture
}

7. Building the Three.js Mesh

Convert the eroded heightfield into a BufferGeometry and compute per-vertex normals from finite differences so the material shades correctly without an extra normal-map pass:

function buildTerrainMesh(field, worldSize) {
  const geo = new THREE.PlaneGeometry(
    worldSize, worldSize, field.size - 1, field.size - 1
  );
  geo.rotateX(-Math.PI / 2);

  const pos = geo.attributes.position;
  for (let i = 0; i < pos.count; i++) {
    const x = i % field.size, z = Math.floor(i / field.size);
    pos.setY(i, field.get(x, z) * 40); // height scale in world units
  }
  geo.computeVertexNormals(); // derives shading normals from the displaced grid

  const mat = new THREE.MeshStandardMaterial({
    vertexColors: true, roughness: 0.95, flatShading: false
  });
  return new THREE.Mesh(geo, mat);
}

Slope-based texturing: use the normal's Y component (dot(normal, up)) as a blend factor between a grass vertex color at flat areas and a rock color on steep slopes — a cheap alternative to a full triplanar shader.

8. Putting It Together

The full generation pipeline runs once, offline (in a Web Worker for a 512² map so the main thread stays responsive), and produces a static mesh:

function generateTerrain(size = 512) {
  const field = new Heightfield(size);

  // 1. Base shape from warped fBm
  for (let z = 0; z < size; z++)
    for (let x = 0; x < size; x++)
      field.set(x, z, warpedHeight(x, z));

  // 2. Alternate hydraulic and thermal erosion in batches
  const erosionParams = {
    maxSteps: 64, inertia: 0.05, minSlope: 0.01,
    sedimentCapacity: 4, depositionRate: 0.3, erosionRate: 0.3,
    gravity: 4, evaporation: 0.02
  };
  for (let batch = 0; batch < 100; batch++) {
    erode(field, 1000, erosionParams);
    thermalErosion(field, 0.6, 0.5);
  }

  // 3. Derive rivers for texturing
  const flow = flowAccumulation(field);

  return { mesh: buildTerrainMesh(field, 1000), field, flow };
}

Performance: 100,000 droplets on a 512² grid take roughly 1-2 seconds on a modern desktop CPU, single-threaded. Move the whole pipeline into a Web Worker and transfer the resulting Float32Array back with a zero-copy Transferable to avoid blocking the render loop.