Simulating Tsunamis: From Shallow-Water Equations to GPU-Accelerated Inundation Maps

Why a wave in a 4,000-metre-deep ocean counts as "shallow water"

Tsunami waves are unusual in that they qualify as shallow-water waves even in the deep ocean, because "shallow" in this context is relative to wavelength, not absolute depth. A tsunami's wavelength can stretch to hundreds of kilometres, vastly larger than even a 4-kilometre-deep ocean basin, so the water column behaves as a thin sheet relative to the wave — justifying models that average flow properties over depth rather than resolving vertical structure. Under this approximation, the wave's phase speed is set by depth alone: c = √(gh), where g is gravitational acceleration and h is water depth. In a typical 4,000 m deep ocean basin this works out to roughly 200 m/s — around 700 km/h, comparable to a jet airliner — which is why a tsunami generated off Japan can reach the west coast of the Americas in under a day.

As the wave approaches the coast and depth h shrinks, c = √(gh) drops accordingly. Since wave energy flux must be roughly conserved, the wave slows and its amplitude grows — the shoaling effect that turns an open-ocean tsunami, often less than a metre high and barely noticeable to a ship, into a multi-metre wall of water by the time it reaches shore.

The governing equations: mass and momentum, depth-averaged

The workhorse model for tsunami propagation is the nonlinear shallow water equations (NSWE), a depth-averaged simplification of the full 3D fluid equations. Working in water depth h(x,y,t), depth-averaged velocity components (u,v), and bathymetry b(x,y), the system conserves mass and momentum:

Mass: ∂h/∂t + ∂(hu)/∂x + ∂(hv)/∂y = 0

Momentum (x-direction): ∂(hu)/∂t + ∂(hu² + ½gh²)/∂x + ∂(huv)/∂y = −gh·∂b/∂x − τₓ

with an equivalent equation in y. The term ½gh² captures hydrostatic pressure, −gh·∂b/∂x is the gravitational forcing from a sloping seabed, and τ represents bottom friction, typically modelled with a Manning-type drag law that grows with velocity and shrinks with depth. These equations are accurate for propagation across the open ocean and are what most operational tsunami forecast models — including NOAA's MOST and the widely used GeoClaw — solve as their core.

Near the coast, where wavelengths shorten and vertical acceleration becomes non-negligible, the shallow-water approximation starts to miss frequency-dependent dispersion effects. Dispersive extensions — Boussinesq-type equations such as the Nwogu (1993) formulation, or the more complete Serre–Green–Naghdi equations — add extra terms that let short, steep waves travel at slightly different speeds than long ones, improving accuracy for nearshore wave shape and for the kind of short-wavelength waves generated by underwater landslides, at meaningfully higher computational cost than the plain shallow-water model.

Numerical methods: capturing a moving bore without blowing up

Solving these equations numerically is complicated by the fact that a tsunami's leading edge can steepen into something resembling a hydraulic bore — a near-discontinuity in water depth — which naive finite-difference schemes handle poorly, producing spurious oscillations or outright instability. Modern solvers instead use finite-volume methods with Godunov-type schemes: the domain is divided into cells, and at each cell interface an approximate Riemann solver computes the flux of mass and momentum given the states on either side. Common choices are the HLL and HLLC solvers, which are robust and computationally cheap, and the Roe solver, which resolves shocks and contact discontinuities more sharply but is more prone to producing unphysical states (like negative water depth) near dry zones — a trade-off between sharpness and robustness that most production solvers resolve in favour of HLLC.

Two further properties matter for stability and correctness. Well-balanced schemes are constructed so that a perfectly still lake sitting over uneven bathymetry stays exactly still numerically — without this so-called C-property (a term introduced by Bermúdez and Vázquez in 1994), the discretisation of the bed-slope source term −gh·∂b/∂x can generate spurious currents purely from numerical error, even with no wave present. Positivity-preserving schemes, combined with high-resolution reconstruction methods like MUSCL or WENO with slope limiters and strong-stability-preserving Runge–Kutta time integration, ensure water depth never goes negative near the wetting/drying front at the coastline — a common failure mode for less careful discretisations. Time steps are limited by the Courant–Friedrichs–Lewy (CFL) condition, typically enforced with a CFL number around 0.4–0.5, based on the fastest wave speed in the domain, |u| + √(gh).

Getting the coastline right: bathymetry, friction, and wetting/drying

A simulation is only as good as the seafloor and coastal elevation data feeding it. Basin-scale propagation can run on relatively coarse bathymetry — grid spacing of a few hundred metres to a few kilometres — but accurate inundation mapping close to shore needs LIDAR-derived topography at metre-scale resolution, since the exact height and shape of the coastline strongly controls how far a wave runs inland. Merging LIDAR survey tiles from different sources or dates is a common source of artificial elevation discontinuities that can visibly distort simulated flood extents if not corrected.

As the wave crosses from wet to dry cells, solvers track a moving wetting/drying front, generally by flagging any cell with depth below a small threshold (on the order of a millimetre to a centimetre) as dry and excluding it from active flux calculations, then reactivating it once inflow raises its depth back above the threshold. Bottom friction is scaled by land-cover type — a Manning roughness coefficient for open grassland is far lower than for dense urban blocks — and buildings or engineered flood barriers are represented either as increased local friction/porosity in coarser models, or as explicit obstacles with overtopping criteria on the finest grids used for site-specific studies.

Making it fast enough to matter: AMR and GPU acceleration

Running an entire ocean basin at metre-scale resolution is computationally infeasible, and mostly unnecessary — most of the deep ocean stays smooth and slowly varying while only the wavefront and the coastal margin need fine detail. Adaptive mesh refinement (AMR), used by solvers such as GeoClaw, addresses this by dynamically refining the computational grid only around the moving wavefront and near coastlines, coarsening it again elsewhere, cutting the total cell count by orders of magnitude compared to a uniform fine grid. Because AMR redistributes work unevenly across the domain, load balancing between compute nodes typically relies on space-filling curves (Morton or Hilbert ordering) to divide the irregular set of active blocks fairly between processors, and conservation across refinement-level boundaries requires careful reflux corrections so mass and momentum aren't lost or created at the interface between coarse and fine regions.

GPU acceleration targets the flux and reconstruction computations that dominate runtime in a finite-volume solver, since these are highly parallel, cell-local operations well suited to GPU architectures. The main practical challenge is keeping data resident on the GPU across timesteps rather than repeatedly transferring it over the relatively slow PCIe bus — GPU-native solvers such as Tsunami-HySEA, built around well-balanced HLLC-type second-order schemes, report large speedups over CPU-only equivalents, which is what makes near-real-time operational forecasting (producing usable inundation estimates within minutes of an earthquake, rather than hours) achievable at all.

Validation: proving the model against reality before trusting it

No tsunami model is used operationally without passing a standard battery of tests. Analytic benchmarks — a perfectly still lake (checking the well-balanced property holds numerically), an idealised dam-break, and a solitary wave running up a sloped beach with a known analytic solution — catch basic implementation bugs. Beyond that, models are checked against the NOAA National Tsunami Hazard Mitigation Program (NTHMP) benchmark problem suite, a standardised set of laboratory and analytic test cases adopted across the tsunami-modelling community specifically so that different research groups' codes can be compared on equal footing. The final and most demanding test is reproducing real historical events — chiefly the 2004 Indian Ocean tsunami and the 2011 Tōhoku tsunami, both extensively instrumented with tide gauges and post-event survey data — comparing simulated arrival times, wave amplitudes, and inundation extents against what was actually recorded. Sensitivity analysis, varying friction coefficients, grid resolution, and earthquake source parameters (often generated from the Okada 1985 analytical model of seafloor displacement from fault rupture), rounds out validation by quantifying how much uncertainty in the inputs translates into uncertainty in the hazard map.