Shallow Water Equations: Build a Tsunami Simulation in an Hour

1. The shallow water assumption

A tsunami in the open Pacific has a wavelength of 100-200 km but the ocean is only ~4 km deep — a depth-to-wavelength ratio of roughly 1:30. Whenever a wave's wavelength is much larger than the fluid depth, the vertical velocity and acceleration inside the fluid column are negligible compared to the horizontal ones. This lets us collapse the full 3D incompressible Navier-Stokes equations (see our companion article on Navier-Stokes in the atmosphere for the general form) into a much simpler 2D system, tracking only water-surface height h(x, y, t) and a depth-averaged horizontal velocity field u(x, y, t) — no vertical dimension at all.

2. The shallow water equations

Dropping vertical structure gives a compact system of three coupled PDEs:

h is the total water column height, u and v are the depth-averaged horizontal velocity components, and g is gravitational acceleration (9.81 m/s²). The right-hand sides are the pressure-gradient force — hydrostatic pressure at any point is proportional to the water column above it, so a sloping surface pushes fluid downhill exactly like a pressure-gradient force in the full Navier-Stokes momentum equation.

For small-amplitude waves, we can drop the non-linear advection terms (u∂u/∂x etc.) entirely, giving the linearised shallow water equations — good enough for open-ocean tsunami propagation, though the full non-linear terms become important as the wave shoals and steepens near the coast.

3. Wave speed and shoaling

Linearising and combining the three equations gives a classic wave equation with speed:

H is the undisturbed water depth. In 4 km deep ocean, c = √(9.81 × 4000) ≈ 198 m/s — over 700 km/h, as fast as a jet airliner, which is why a tsunami generated by an earthquake can reach a coastline thousands of kilometres away within hours. Because c depends on depth, as the wave approaches shore and H decreases, the wave slows down.

Conservation of energy flux (roughly, energy flux ∝ amplitude² × group velocity must stay constant) means that as speed drops, amplitude must rise to compensate — a phenomenon called shoaling. Green's law gives a common approximation:

A tsunami that's a barely-noticeable 0.5 m swell in 4 km deep ocean can grow to several metres by the time H drops to a few tens of metres near shore — exactly what the shallow water equations predict, and exactly why an earthquake felt at sea produces a devastating wave on land.

4. Finite-difference discretisation

To simulate this on a computer, we discretise h, u and v onto a regular 2D grid with spacing dx, and replace spatial derivatives with central differences:

and advance the grid forward in time with a simple explicit (forward Euler) time-step: compute all spatial derivatives at the current time, use them to compute the rate of change of h, u, v, then step everything forward by dt.

5. The CFL stability condition

Explicit schemes are only stable if information can't cross more than one grid cell per time-step. Since the fastest signal in this system travels at the wave speed c = √(gH), the Courant-Friedrichs-Lewy (CFL) condition requires:

Pick dt any larger and the simulation will diverge into NaN within a few frames — this is the single most common bug when first implementing a shallow-water solver.

6. Implementation in JavaScript

class ShallowWaterSim {
  constructor(nx, ny, dx) {
    this.nx = nx; this.ny = ny; this.dx = dx;
    this.g = 9.81;
    this.H  = new Float32Array(nx * ny).fill(50);  // still-water depth (bathymetry)
    this.h  = new Float32Array(nx * ny);              // surface elevation (perturbation)
    this.u  = new Float32Array(nx * ny);
    this.v  = new Float32Array(nx * ny);
  }

  idx(i, j) { return j * this.nx + i; }

  // Drop a Gaussian bump — simulates an earthquake-generated initial displacement
  addDisturbance(cx, cy, amplitude, radius) {
    for (let j = 0; j < this.ny; j++)
      for (let i = 0; i < this.nx; i++) {
        const d = Math.hypot(i - cx, j - cy);
        this.h[this.idx(i, j)] += amplitude * Math.exp(-(d * d) / (2 * radius * radius));
      }
  }

  step(dt) {
    const { nx, ny, dx, g, H, h, u, v } = this;
    const hNew = new Float32Array(h.length);
    const uNew = new Float32Array(u.length);
    const vNew = new Float32Array(v.length);

    for (let j = 1; j < ny - 1; j++)
      for (let i = 1; i < nx - 1; i++) {
        const k = this.idx(i, j);
        const total = H[k] + h[k]; // total column height

        // Central differences for gradients
        const dhdx = (h[this.idx(i+1,j)] - h[this.idx(i-1,j)]) / (2 * dx);
        const dhdy = (h[this.idx(i,j+1)] - h[this.idx(i,j-1)]) / (2 * dx);

        // Momentum: du/dt = -g * dh/dx  (linearised, no advection)
        uNew[k] = u[k] - g * dhdx * dt;
        vNew[k] = v[k] - g * dhdy * dt;

        // Continuity: dh/dt = -div(total * u)
        const fluxRight = total * u[this.idx(i+1,j)];
        const fluxLeft  = total * u[this.idx(i-1,j)];
        const fluxUp    = total * v[this.idx(i,j+1)];
        const fluxDown  = total * v[this.idx(i,j-1)];
        const div = (fluxRight - fluxLeft) / (2 * dx) + (fluxUp - fluxDown) / (2 * dx);
        hNew[k] = h[k] - div * dt;
      }

    this.h = hNew; this.u = uNew; this.v = vNew;
    this.applyBoundaries();
  }
}

7. Reflecting boundaries and a coastline

A domain edge needs a boundary condition or the finite-difference stencil reads out of bounds. The simplest physically reasonable choice is a reflecting (solid wall) boundary: zero the normal velocity component and mirror the height gradient, so waves bounce back instead of vanishing or wrapping around.

applyBoundaries() {
  const { nx, ny, u, v, h } = this;
  for (let j = 0; j < ny; j++) {
    u[this.idx(0, j)]      = 0; u[this.idx(nx - 1, j)] = 0;
    h[this.idx(0, j)]      = h[this.idx(1, j)];
    h[this.idx(nx - 1, j)] = h[this.idx(nx - 2, j)];
  }
  // (repeat symmetrically for the top/bottom edges using v)
}

To simulate a beach and watch shoaling in action, make the still-water depth H shallower toward one edge of the grid instead of constant — a linear ramp from 50 m down to 0.5 m over the last 20% of the domain is enough to see the wave slow down and visibly grow in amplitude exactly as Green's law predicts, right before it "breaks" as the linear approximation breaks down at very shallow depth.