Shallow Water Equations: How Tsunamis Cross an Ocean

1. Why "Shallow" Works for an Ocean-Spanning Wave

"Shallow" sounds like an odd word to describe the Pacific Ocean, but it refers to a ratio, not an absolute depth. A wave is "shallow-water" whenever its wavelength λ is much larger than the water depth h — formally, when h/λ ≪ 1 (in practice, roughly h/λ < 1/20).

A tsunami generated by an undersea earthquake typically has a wavelength of 100–500 km, while the average ocean depth is only about 4 km. The ratio h/λ is therefore around 0.01–0.04 — comfortably in shallow-water territory, even though the water itself is several kilometres deep. Wind-driven surface waves, by contrast, have wavelengths of tens to hundreds of metres and are almost always deep-water waves, governed by completely different physics.

This distinction matters because in the shallow-water regime, the horizontal length scale so dominates the vertical one that the vertical momentum equation collapses to simple hydrostatic balance, and the flow can be treated as essentially two-dimensional — a huge simplification that makes basin-scale tsunami simulation computationally tractable.

2. Deriving the Shallow Water Equations

The shallow water equations (SWE) are obtained by depth-averaging the incompressible Navier-Stokes equations under the hydrostatic assumption. Let h(x,y,t) be the water depth, η(x,y,t) the free-surface elevation above a reference level, b(x,y) the bed elevation (so h = η − b), and u = (u,v) the depth-averaged horizontal velocity. Integrating mass and momentum conservation from the bed to the free surface gives:

Expanded into components, with g the gravitational acceleration:

The term ½gh² arises from integrating the hydrostatic pressure P = ρg(η − z) through the water column — it plays exactly the role that pressure P plays in the full Navier-Stokes momentum equation, except now it acts like a nonlinear "pressure" that depends on depth squared. This is why the SWE are also called the Saint-Venant equations (introduced by Adhémar Barré de Saint-Venant in 1871 for open-channel river flow, decades before they were applied to tsunamis).

The right-hand-side terms represent: −gh∇b, the force from a sloping seabed (bathymetry); fv h and −fu h, the Coriolis force from Earth's rotation (f = 2Ω sin φ, important for basin-scale, multi-hour propagation); and −c_f u|u|, quadratic bottom friction, negligible in the deep ocean but critical near the coast.

3. Wave Speed: c = √(gh)

Linearising the shallow water equations around a still ocean of depth h (small disturbance η ≪ h, no advection or friction terms) yields the classic linear wave equation for the surface elevation:

which is satisfied by waves travelling at the shallow-water phase speed:

This single formula explains why a tsunami is one of the fastest natural phenomena on the planet. In the open Pacific, where h ≈ 4,000 m:

That is comparable to a commercial jet's cruising speed — a tsunami can cross the entire Pacific Ocean in less than a day. Crucially, in shallow-water theory c depends only on depth, not on wavelength — the wave is non-dispersive: every Fourier component of the initial disturbance travels at the same speed, so the wave keeps its shape as it propagates across the open ocean instead of spreading out.

4. How an Earthquake Generates a Tsunami

Roughly 80% of tsunamis are triggered by subduction-zone earthquakes. When one tectonic plate suddenly slips beneath another, a patch of seafloor — often hundreds of kilometres long — is vertically displaced by metres in a matter of seconds. Because water is nearly incompressible and the shallow-water assumption treats the vertical velocity as instantaneously transmitted to the surface, the initial condition for a tsunami simulation is simple:

where Δb is the vertical seafloor deformation computed from an elastic dislocation model of the fault rupture (the widely used Okada (1985) formulas give surface deformation from fault geometry, slip, and moment magnitude). This seafloor "bump" becomes the initial free-surface elevation, and the SWE are then integrated forward in time to propagate it.

Other tsunami sources fit the same initial-value framework with a different Δb or initial velocity field: submarine landslides (a rapidly moving mass redistributes water — sometimes with much shorter wavelengths, invalidating pure shallow-water assumptions), volcanic flank collapses, and — rarely — asteroid impacts or explosive volcanic eruptions like Hunga Tonga in 2022, which generated atmospheric pressure waves that also coupled into the ocean.

5. Shoaling, Green's Law, and Run-up

As a tsunami approaches the coast, depth h decreases, and c = √(gh) decreases with it. Since wave energy flux must be approximately conserved along a slowly varying channel or shelf, the wave must grow in height to compensate for its slowing speed. This amplification process is called shoaling.

For a wave travelling into gradually decreasing depth (no reflection, energy flux conserved), Green's law gives the height scaling:

Take a tsunami of height H₁ = 0.5 m in h₁ = 4,000 m of water shoaling onto a coastal shelf of h₂ = 20 m:

Green's law alone already gives a factor of ~4 amplification, and it breaks down entirely in the final metres of run-up, where nonlinear effects, bottom friction, and coastline geometry (bays, river mouths, and headlands can funnel and further amplify the wave) dominate. Real run-up heights routinely exceed the open-shelf estimate, which is why some 2011 Tōhoku tsunami run-up measurements reached nearly 40 m in narrow inlets, even though the offshore wave height was on the order of a few metres.

The Leading-Edge Trough

Many tsunamis arrive coastal-first as a trough rather than a crest — the sea visibly recedes, exposing the seabed, before the main crest arrives minutes later. This happens whenever the initial seafloor deformation has a leading depression adjacent to the main uplift, or through wave dispersion; it is one of nature's few reliable, if terrifying, warning signs and is the reason coastal evacuation guidance says: if the sea drastically recedes, do not go look — run to high ground immediately.

6. Dispersion and the Limits of the Shallow Water Model

The classical shallow water equations are derived by assuming the pressure is perfectly hydrostatic, which is equivalent to dropping all vertical acceleration. This is an excellent approximation for h/λ ≪ 1, but becomes inaccurate for shorter-wavelength disturbances such as landslide-generated tsunamis, meteotsunamis, or the leading edge of an earthquake tsunami as it steepens near shore.

A more complete linear theory (Airy wave theory) gives the exact dispersion relation for water waves of any depth:

where ω is angular frequency and k = 2π/λ is wavenumber. Two useful limits:

• Shallow water (kh ≪ 1): tanh(kh) ≈ kh, giving ω² ≈ g k² h, i.e. phase speed c = ω/k → √(gh) — exactly the non-dispersive shallow-water result.
• Deep water (kh ≫ 1): tanh(kh) → 1, giving c = √(g/k) — speed now depends on wavelength, and waves disperse (this is the regime of ordinary wind-driven ocean swell).

To capture weak dispersion while retaining computational efficiency, tsunami modellers use Boussinesq-type equations, which add higher-order derivative correction terms to the SWE to approximate the tanh(kh) behaviour without solving the full 3D free-surface Navier-Stokes problem. These are essential for accurately modelling landslide tsunamis and the fine structure of wave trains after long-distance propagation.

7. Numerical Methods: Finite Volumes and Well-Balanced Schemes

The SWE are a system of nonlinear hyperbolic conservation laws, closely related in structure to the compressible Euler equations of gas dynamics (depth h plays the role of density, and ½gh² plays the role of pressure). This means the same numerical machinery developed for shock-capturing in compressible flow — Godunov-type finite volume methods — applies directly to tsunami simulation, and bores/hydraulic jumps in shallow water are the exact analogue of shock waves in gas dynamics.

A typical finite volume update on a computational cell i, over a discrete timestep Δt, takes the conservative form:

The numerical flux F at each cell interface is typically computed with an approximate Riemann solver — HLL (Harten-Lax-van Leer) or Roe's scheme are the most common choices, chosen for their ability to resolve the bores that form at a tsunami's leading edge as it steepens near shore.

Well-Balanced Schemes and Wetting-Drying

A naive discretisation of the bathymetry source term −gh∇b fails a critical sanity check: a perfectly still ocean over variable bathymetry should stay perfectly still. Schemes that preserve this "lake-at-rest" steady state exactly are called well-balanced schemes, and they are mandatory for tsunami codes — otherwise numerical noise from the bathymetry gradient alone can spuriously generate spurious waves large enough to mask the real signal.

Coastal inundation additionally requires wetting-and-drying algorithms: cells where h → 0 must be handled carefully (typically with a thin-film threshold and special flux limiting) so the simulation can represent the shoreline advancing inland during run-up and retreating during draw-down without the depth going negative or the velocity blowing up.

8. Real Tsunami Warning Systems

Operational tsunami warning centres — such as the US NOAA Pacific Tsunami Warning Center and Japan's JMA — run shallow water models built on exactly this framework, but under extreme time pressure: a forecast must be issued within minutes of an earthquake, long before the tsunami itself arrives at any coast.

The standard operational pipeline looks like this:

• Seismic detection — global seismometer networks estimate earthquake location, depth, and moment magnitude within 1–3 minutes.
• Source modelling — an Okada-type fault model converts the earthquake parameters into an initial seafloor deformation Δb.
• Precomputed scenario database — because full SWE simulation is too slow to run from scratch in real time, centres like NOAA maintain databases of thousands of precomputed unit-source tsunami propagations (the SIFT / MOST system) that are linearly combined to match the actual earthquake source in seconds.
• DART buoy confirmation — deep-ocean bottom-pressure sensors (Deep-ocean Assessment and Reporting of Tsunamis buoys) measure the tiny pressure signature of the passing wave in open water and are used to correct and validate the forecast in near-real time.
• Coastal forecast and inundation mapping — high-resolution, well-balanced SWE runs with detailed local bathymetry then estimate arrival time, wave height, and inundation extent for each populated coastline segment.

The same physics that makes GPU-accelerated real-time water simulations possible for games and visualisation — a depth-averaged grid updated every frame with a finite-volume or finite-difference shallow water solver — is, at much larger scale and with life-safety stakes, the exact computation running inside a tsunami warning centre the moment a subduction earthquake is detected.