Shallow Water Equations: 2D Height & Velocity Field Solver

The 1D version tracks a single cross-section h(x,t) and u(x,t) along one line. This 2D version solves the full coupled system over a grid — h(x,y,t), u(x,y,t) and v(x,y,t) — so waves can spread outward, curve, and interfere the way real water waves do, which a 1D solver cannot represent.

Wave speed in shallow water is c = √(gh): taller, deeper regions of the wave travel faster than shorter, shallower regions ahead of them. As a wave's front steepens enough, the local surface gradient crosses a threshold and the cell is counted as a "breaking-front" cell — a numerical stand-in for a bore or a breaking wave.

Its bed depth decreases from left to right. As the initial pulse moves into shallower water, its speed c = √(gh) drops, and because momentum flux must still be conserved, the amplitude grows to compensate — the classic shoaling amplification that makes ocean swells rear up and break as they approach a beach.

Two fixed grid cells are driven with an oscillating height, each radiating circular waves like ripples from two stones dropped together. Where the two wave trains arrive in phase they add constructively into taller crests; where they arrive out of phase they cancel — the same principle behind double-slit interference, but for real nonlinear water waves.

The solver is written in conservative flux form, so the total water volume Σh·dx·dy over the whole grid should stay close to its initial value for a closed, wall-bounded domain like this one. Tracking it is a running check that the finite-difference scheme isn't numerically creating or destroying water.