Tsunami propagation (finite-volume SWE, HLL + minmod)

Platform v3.0.0 Docs v3.0.0 MEASURED LOCALLY

Full description: Read the full narrative

A tsunami scenario starts from a finite-fault seafloor displacement (simplified Okada-1985-in-spirit rectangular dip-slip patch with moment-magnitude scaling) and propagates the free surface with a nonlinear shallow-water finite-volume solver (MUSCL/minmod reconstruction, HLL fluxes, semi-implicit Manning friction) over real GEBCO 2020 bathymetry sampled server-side; a synthetic bathymetry fallback is explicitly disabled.

Execution mode: Kaggle kernel (numpy-only, CPU accelerator instance) server/kaggle/simRunner.ts:126-141 (not in LOCAL_SIM_TYPES) kaggle-kernels/tsunami-sim/main.py:228-229
Verification: SOLVER MEASURED locally in-process (conservation gate PASS; flat-bed propagation stable; REAL-GEBCO runs diverge — see Known limitation)
Wire scenario type: tsunami_wave — contract at src/services/kaggleSim.ts:28-240; job API at Simulation job API

Evidence legend

MEASURED marks values produced by executing the actual code on this machine (commands in Reproduction; raw logs summarised in Methodology → Evidence runs).

Governing equations and coefficients

Source attribution as stated in the kernel header (verbatim):

Initial source: simplified Okada-style elastic dislocation uplift — vertical seafloor displacement from finite-fault slip transferred to the water surface (Okada 1985 in spirit; here a rectangular fault patch with dip-slip lobes, and moment-magnitude scaling).
kaggle-kernels/tsunami-sim/main.py:10-13

∂η/∂t + ∇·(H·u) = 0; ∂u/∂t + … = −g·∇η; c = √(g·H)

Header form (H = still-water depth + η). The implemented solver is the NONLINEAR conservative SWE (h·u² and g·h²/2 fluxes), not the linearised form. kaggle-kernels/tsunami-sim/main.py:7-9 kaggle-kernels/tsunami-sim/swe_solver.py:64-75

M0 = 10^(1.5·M + 9.1) N·m; amp = seafloor_displacement_m · √(M0/M0_ref) / 5.6, M0_ref = 10^(1.5·8+9.1)

Moment scaling of the surface displacement (normalised to the requested displacement at M 8). kaggle-kernels/tsunami-sim/main.py:262-267

L_half = 20·10^(0.5(M−8)) km; W_half = max(0.4·L_half, 8) km; slip = clip(1−(along/L)²,0,1)·clip(1−(across/W)²,0,1)·tanh(across/(0.5·W))

Rupture dimensions ("Wells & Coppersmith style") and dip-slip uplift/subsidence lobes; strike default 135°; applied to ocean cells (bathy > 50 m) only. kaggle-kernels/tsunami-sim/main.py:274-287,260,94

Face flux: MUSCL minmod reconstruction; HLL with sl = min(u−c), sr = max(u+c); dt = 0.9·dx/(√2·c_max)

Update scheme and main-loop timestep (the solver's own estimate_cfl_dt uses safety 0.45 — the kernel main path passes its own dt). kaggle-kernels/tsunami-sim/swe_solver.py:34-36,77-113

friction = g·n²·|u|/H^(4/3), clipped ≤ 0.05 s⁻¹, applied as (1+dt·friction)⁻¹; n = 0.03

Semi-implicit Manning damping. kaggle-kernels/tsunami-sim/swe_solver.py:200-206

∂(hu)/∂t += −dt·g·H·∂b/∂x (central difference)

Explicit bathymetry source term — NOT hydrostatically re-balanced per face; see Known limitation. kaggle-kernels/tsunami-sim/swe_solver.py:188-193

Coefficients and named constants

Constants and defaults.
Name	Value	Meaning	Source
g	9.81 m/s²	Gravity	kaggle-kernels/tsunami-sim/main.py:232
manning_n	0.03	Basal roughness	kaggle-kernels/tsunami-sim/main.py:307 kaggle-kernels/tsunami-sim/swe_solver.py:163
ocean threshold	bathy > 50 m	Source applied only on ocean cells	kaggle-kernels/tsunami-sim/main.py:94,254
default dx	5 km (when no extent_km)	Cell size	kaggle-kernels/tsunami-sim/main.py:233-234
dt (default)	0.9·dx/(√2·√(g·H_max))	CFL-style timestep from deep-water celerity	kaggle-kernels/tsunami-sim/main.py:236-238
strike default	135°	Rupture orientation	kaggle-kernels/tsunami-sim/main.py:260

Parameter contract

Tsunami parameters. Shared request envelope (validated at the client boundary before dispatch).

Tsunami parameters.
UI control	Wire field	Unit	Valid range (UI · wire · kernel)	Default	Physical meaning	Source
Earthquake Magnitude slider	`magnitude`	M	UI 5–9.5 step 0.1; wire 5–9.7	UI 8.5 · kernel 8.5	Sets M0, rupture L/W and amplitude	src/components/scenarios/ScenarioEditor.tsx:109-113 src/services/kaggleSim.ts:113 kaggle-kernels/tsunami-sim/main.py:191,262-276
Max Wave Height slider	`seafloor_displacement_m`	m	UI 1–30 m (hidden derived value preferred: clamp(1, 40, 5·10^(M−6))); wire >0, ≤ 40	derived at M 8.5 → 40 m; kernel 5	Near-field vertical displacement reference amplitude	src/components/scenarios/ScenarioEditor.tsx:109-113 src/services/kaggleSim.ts:114,473,660 kaggle-kernels/tsunami-sim/main.py:192
Hypocenter Depth slider	— (not sent)	km	UI 5–100	20	Panel context + arrival-time hint only; long-wave physics is depth-independent	src/components/scenarios/ScenarioEditor.tsx:109-113 src/services/kaggleSim.ts:663-664
duration_minutes	`duration_minutes`	min	wire ≤ 720; builder fixes 30	kernel 30	Propagation window	src/services/kaggleSim.ts:115,474 kaggle-kernels/tsunami-sim/main.py:193
Bathymetry	`bathymetry(+_gs)` → `bathy_b64/min/span`	m positive-down	Required: run raises without real GEBCO	auto server-sampled	GEBCO 2020 via api.opentopodata.org/v1/gebco2020 (64² default)	src/services/kaggleSim.ts:116-124 server/kaggle/bathymetry.ts:10,28 kaggle-kernels/tsunami-sim/main.py:227-230
benchmark_stations (optional)	array ≤16	—	row/col + observed eta(t) series	—	Compare modeled η(t) with observed gauge series	src/services/kaggleSim.ts:125-133

Outputs

Files written per run (also served by `GET /api/kaggle/simulate/:jobId/grid/:name` and `…/geotiff/:name`)
Field	Unit	Meaning	Source
`water_height.npy`	m (η)	Final free-surface elevation grid	kaggle-kernels/tsunami-sim/main.py:421-425
`bathymetry.npy`	m	Depth grid used (positive-down)	kaggle-kernels/tsunami-sim/main.py:421-425
`snapshots_height.npy` + `snapshot_times.npy`	m; min	Propagation frames	kaggle-kernels/tsunami-sim/main.py:421-425
`metadata.json` + `validation.json` + `benchmark.json`	—	model "okada_finite_fault_dislocation", solver "finite_volume_swe_hll_minmod", final stats, benchmark comparison	kaggle-kernels/tsunami-sim/main.py:336-364,430-438

Spatial-origin semantics

Square domain decoded from the uint16-quantized GEBCO sample (bathy_min + u16/65534·bathy_span, sentinel 65535 = NaN→min) and bilinearly upsampled to the run grid; positive-down depth, negative = dry land; row 0 = north, matching the server-side sampler that mirrors the client's cell-centre geometry. The rupture is placed in the deep-ocean region of the box, not the centre. kaggle-kernels/tsunami-sim/main.py:208-226,59,79-99 server/kaggle/bathymetry.ts:18-24

Documented validity limits (verbatim from the code)

'Real GEBCO bathymetry is required for tsunami simulations; synthetic fallback is disabled.' · the linear-in-form header describes c = sqrt(g·H) propagation; implementation integrates the nonlinear conservative SWE · internal conservation proof states: 'Flat bathymetry (constant depth) — no source terms, pure wave propagation.'
kaggle-kernels/tsunami-sim/main.py:228-229 kaggle-kernels/tsunami-sim/main.py:7-9 kaggle-kernels/tsunami-sim/swe_solver.py:235

main() runs the closed-box conservation proof (flat 2000 m bathy, 500 steps, tol 1e-9) before the simulation and refuses non-conservative runs: "Tsunami solver failed closed-box conservation check (|Δ|=…). Refusing to run a non-conservative simulation." kaggle-kernels/tsunami-sim/main.py:406-417

Measured runs and timings

Executed on an Apple-Silicon laptop (Python 3.14.5, NumPy 2.4.2). These are observations of one environment, not performance guarantees.
Run / check	Measured result
Conservation gate (re-executed)	[VERIFY] closed-box gs=64 steps=500: init=8194010.359005 final=8194010.359005 \|Δ\|=0.00e+00 → PASS
Flat-bed propagation (4000 m constant bathy, 96²/300 km, dt = 3 s → CFL 0.27, 20 min)	Stable for the duration: 1.21 m uplift → max η 0.64 m radiated; energy 0.07 GJ; wall 0.6 s
Flat bed, kernel default dt (96²/300 km, dt = 0.9·dx/(√2·c) = 10.04 s, printed “CFL 0.90 (OK)”, 20 min)	Diverges (max \|η\| ~ 10²³⁰ m): the kernel’s own dt formula exceeds the ~0.5 stability bound of second-order (MUSCL/HLL) updates; its estimate_cfl_dt uses safety 0.45 but main() passes its own dt
Flat bed, small box (48² at dt = 3 s)	Stable for ~40 steps, then monotonic growth from periodic wrap-around coupling at the roll-based face divergence: 0.27 → 22 m by 20 min
Sloped synthetic bed (dt = 3 s)	Diverges at low CFL — the explicit bathymetry source drives unbounded growth on any slope
REAL GEBCO 2020 (16×16 sampled off Tohoku, 38.5°N 143.5°E, via the same api.opentopodata.org/v1/gebco2020 endpoint the server uses; M 8.5, 5 m, 20 min, 300 km box)	Diverges at both default dt (7.25 s, “CFL 0.90 (OK)”) and dt = 3 s (CFL 0.37): final max_wave ~10³⁰⁰ m, E = nan

Known numerical limitation

Tsunami solver — non-flat bathymetry instability

This review reproduced three distinct unbounded-growth mechanisms in the standalone kernel: (1) the default main-loop timestep uses a 0.9 CFL that second-order MUSCL/HLL updates cannot satisfy (flat bed still diverged at CFL 0.90; the solver's own estimate_cfl_dt applies safety 0.45 but main() overrides it); (2) the roll-based face divergence wraps opposite domain edges together, so waves grow after boundary interaction even on a flat bed in a small box; (3) on any non-flat bathymetry — including a real GEBCO 2020 sample — the explicit central-difference source term drives divergence at every tested timestep (max |η| → 10^300 m, energy NaN). The closed-box conservation proof deliberately uses a flat bed where the source vanishes, so it cannot detect (2) or (3). Until the scheme gains a well-balanced source treatment (e.g. the Audusse reconstruction the flood kernel uses) and a transmissive boundary, tsunami output must be treated as UNVALIDATED. Production runs go through Kaggle where this review could not execute them.

Reproduction

Commands

# Conservation gate + flat-bed propagation + diverging real-bed case:
cd /tmp && mkdir ts && cp <repo>/kaggle-kernels/tsunami-sim/{main.py,swe_solver.py} ts/ && cd ts
python3 - <<'PY'
import base64, json, numpy as np, subprocess, os
v = np.full((16,16), 4000.0)          # flat bed  — stable
# v = (np.mgrid[0:16,0:16][1]/15)*-4000+4000   # sloped bed — diverges
span=v.max()-v.min(); u16=np.round((v-v.min())/max(span,1e-9)*65534).astype('<u2')
json.dump({'grid_size':96,'magnitude':8.5,'seafloor_displacement_m':5,'duration_minutes':20,
  'extent_km':300,'dt_s':3.0,'bathy_gs':16,'bathy_b64':base64.b64encode(u16.tobytes()).decode(),
  'bathy_min':float(v.min()),'bathy_span':float(span)}, open('params.json','w'))
PY
TERRANOETIS_OUT_DIR=$PWD python3 main.py | grep -E "VERIFY|DONE"   # writes fail off-Kaggle: hard-coded /kaggle/working (main.py:421)