[Paper Review] Beating the Storm Clock — Casulli–Walters (2000) Semi-Implicit Shallow Water on Unstructured Orthogonal Grids

Celerity leaves the stability bound only while the free-surface matrix keeps a diagonal surplus of one.

In June 1990, Hamburg's BAW lab had to replay 13 days of tide in the Jade–Weser. 150 thousand triangles. A 5-minute step. Faster than real time, or the next forecast never starts.

An explicit step lets free-surface celerity $c=\sqrt{gh}$ set the CFL bound. At 28 m, $c \approx 16.6$ m/s. A shortest edge near $\sqrt{10}$ m leaves 0.2 s. Thirteen days is 1,123,200 s, so 0.2 s means 5,616,000 steps — the paper's 5.6 million. Five minutes means 3,744 steps. The ratio is exactly 1,500, because 300 / 0.2 does not care which fortnight you picked.

Casulli and Walters (2000) keep that ratio alive by leaving only the free-surface gradient at the new time. Advection stays at the old time.

Drag the CFL slider past 1 and the red explicit curve blows up, while raising θ only shrinks the cyan amplitude.

CFL ratio (Δt / Δt_explicit): 2.00θ (implicitness): 0.60

CFL > 1 ⇒ explicit blows up; semi-implicit stays bounded for any CFL when θ ≥ 1/2.

At CFL=6 the cyan crest is still there. Push θ to 1 and that crest dies quickly; the band where amplitude survives in practice is 0.55 to 0.6.

Paper metadata#

Authors: Vincenzo Casulli (Univ. of Trento), Roy A. Walters (USGS)
Journal: International Journal for Numerical Methods in Fluids, 32, 331–348 (2000)
Title: "An unstructured grid, three-dimensional model based on the shallow water equations"
Keywords: semi-implicit, shallow water, unstructured orthogonal grid, wetting/drying

Why an explicit solver fails on estuaries#

At 10 m, $c \approx 10$ m/s. At 100 m, $c \approx 31$ m/s. Flow speed $u$ usually sits under 1 m/s. The tide crosses a reach in $L/u$ , but the explicit step is stuck at $\Delta x/c$ . With $c/u \approx 10\sim30$ , the code slices time that much finer than the water needs.

An unstructured mesh dumps the whole limit onto the shortest edge. Big Lost River is 1.2 km by 1 km, 12,622 triangles, shortest side near 5 m, explicit limit about 0.5 s. Divide 5 m by 0.5 s and you get 10 m/s. The deep channel's 16 m/s is not what stops the clock. The short edge in shallow water does.

Slow terms stay explicit. The fast one goes implicit. Once only the free-surface pressure gradient is implicit, $\sqrt{gh}$ leaves the stability bound. What remains is a limit from $u$ .

Orthogonal unstructured grids — the Voronoi–Delaunay pair#

Plain finite differences on an unstructured mesh break in two places. Cell spacing changes with shape, so the gradient leans. A curvilinear escape adds transformation terms that cut accuracy and stability together.

Casulli and Walters restricted the mesh instead. An orthogonal unstructured grid: the segment between neighboring centers is perpendicular to the shared edge. A Voronoi polygon and its dual Delaunay triangulation meet that test when every triangle is acute. One obtuse angle and the circumcenter leaves the triangle. The link of centers can stay perpendicular to the edge and still miss the cell interior. How long point insertion keeps the circumcircle empty is how long this difference stays valid.

Orthogonality buys one difference. With only the edge-normal velocity $u_j$ , $\partial \eta / \partial n$ is $(\eta_{i_2} - \eta_{i_1})/\delta_j$ . No curvature fix.

Semi-implicit splitting: the θ method untangles two time scales#

The θ method splits time with an implicitness factor $\theta \in [1/2, 1]$ . Momentum sees the free-surface gradient as $\theta \eta^{n+1} + (1-\theta)\eta^n$ . The free-surface equation sees velocity as $\theta u^{n+1} + (1-\theta)u^n$ .

Discretized momentum, vertical viscosity and wind stress included:

\mathbf{A}_j^n \mathbf{U}_j^{n+1} = \mathbf{G}_j^n - \theta g \frac{\Delta t}{\delta_j}\!\left[\eta_{i(j,2)}^{n+1} - \eta_{i(j,1)}^{n+1}\right]\Delta\mathbf{Z}_j^n

$\mathbf{U}_j^{n+1}$ is the normal velocity on edge $j$ . $\mathbf{A}_j^n$ is the tridiagonal of vertical viscosity, wind, and bottom friction. $\mathbf{G}_j^n$ holds advection, Coriolis, and horizontal friction. Each edge is its own $N_z \times N_z$ system, and there are $N_s$ of them.

The free-surface equation:

P_i \eta_i^{n+1} = P_i \eta_i^n - \theta \Delta t \sum_{l=1}^{S_i} s_{i,l}\lambda_{j(i,l)} [\Delta\mathbf{Z}^\top \mathbf{U}]_{j(i,l)}^{n+1} - (1-\theta)\Delta t \sum_{l=1}^{S_i} s_{i,l}\lambda_{j(i,l)} [\Delta\mathbf{Z}^\top \mathbf{U}]_{j(i,l)}^{n}

$P_i$ is the cell area, $\lambda_j$ the edge length, $s_{i,l}$ the outward sign.

Wave-equation reduction and the SPD payoff#

Left coupled, the pair has $N_z N_s + N_p$ unknowns. Solve the momentum row for $\mathbf{U}_j^{n+1}$ as a function of $\eta^{n+1}$ and substitute. If $\mathbf{A}$ is positive definite, so is its inverse.

What drops out is a discrete wave equation.

P_i \eta_i^{n+1} - g \Delta t^2 \theta^2 \sum_{l=1}^{S_i} \frac{s_{i,l}\lambda_{j(i,l)}}{\delta_{j(i,l)}}\!\left[(\Delta\mathbf{Z})^\top \mathbf{A}^{-1}\Delta\mathbf{Z}\right]_{j(i,l)}^{n}\!\!\left(\eta^{n+1}_{i[j(i,l),2]} - \eta^{n+1}_{i[j(i,l),1]}\right) = R_i

One unknown, $\eta^{n+1}$ . Size $N_p$ . Symmetric, strongly diagonally dominant, positive definite. PCG fits. The paper reports convergence in tens of iterations even on large meshes.

On the 1D toy there is no vertical structure, so $\mathbf{A}$ is 1 and the code's $k = g H \Delta t^2 \theta^2 / \Delta x^2$ is the neighbor coefficient. At four times the explicit limit, $c^2 = gH$ forces $k = 16\theta^2$ . At θ=0.6, $k=5.76$ . The interior diagonal $1+2k$ is 12.52. The two off-diagonal entries sum to 11.52. The surplus is always 1, and that 1 is the storage term $P_i \eta$ . Grow the step and $k$ sticks to the diagonal and the off-diagonal together, so the surplus does not grow. At the slider's CFL=6 and θ=0.6, $k=12.96$ , the diagonal is 26.92, and the off-diagonal sum is 25.92. Surplus: still 1. A more implicit step does not make the matrix safer. What dies is the amplitude.

Once $\eta$ is known, $\mathbf{U}_j^{n+1}$ falls back to a small tridiagonal on each edge. Vertical velocity is recovered from discrete continuity.

A 1D version in Python#

Same split. $N=200$ , $L=1$ , $H=1$ , $g=9.81$ , θ=0.6.

import numpy as np
from scipy.linalg import solve_banded
 
class CasulliShallowWater1D:
    """1D linear shallow water — semi-implicit θ method."""
 
    def __init__(self, N=200, L=1.0, H=1.0, g=9.81, theta=0.6):
        self.N, self.L, self.H, self.g, self.theta = N, L, H, g, theta
        self.dx = L / N
        self.eta = np.zeros(N)            # free-surface at cell centers
        self.u = np.zeros(N + 1)          # face-normal velocity, walls = 0
 
    def initial_bump(self, x0=0.3, sigma=0.05, amp=0.05):
        x = (np.arange(self.N) + 0.5) * self.dx
        self.eta = amp * np.exp(-((x - x0) ** 2) / (2 * sigma ** 2))
        self.u[:] = 0.0
 
    def step(self, dt):
        N, dx, g, H, th = self.N, self.dx, self.g, self.H, self.theta
        k = g * H * dt * dt * th * th / (dx * dx)   # wave-equation coefficient
        alpha = H * dt / dx
 
        # G_j: explicit part collected at edges (walls = 0)
        G = np.zeros(N + 1)
        G[1:N] = self.u[1:N] - g * (dt / dx) * (1 - th) * (self.eta[1:] - self.eta[:-1])
 
        # Tridiagonal: -k η_{i-1} + (1+2k) η_i - k η_{i+1} = R_i
        ab = np.zeros((3, N))
        ab[0, 1:] = -k                         # super-diagonal
        ab[2, :-1] = -k                        # sub-diagonal
        ab[1, :] = 1.0 + 2 * k
        ab[1, 0] = 1.0 + k                     # left wall
        ab[1, -1] = 1.0 + k                    # right wall
 
        rhs = self.eta.copy()
        rhs -= alpha * th * (G[1:] - G[:-1])
        rhs -= alpha * (1 - th) * (self.u[1:] - self.u[:-1])
 
        new_eta = solve_banded((1, 1), ab, rhs)
 
        new_u = self.u.copy()
        new_u[1:N] = self.u[1:N] - g * (dt / dx) * (
            th * (new_eta[1:] - new_eta[:-1]) + (1 - th) * (self.eta[1:] - self.eta[:-1])
        )
        new_u[0] = new_u[-1] = 0.0
        self.eta, self.u = new_eta, new_u
        return self.eta
 
if __name__ == "__main__":
    sim = CasulliShallowWater1D(N=200, theta=0.6)
    sim.initial_bump()
    c = np.sqrt(sim.g * sim.H)
    dt_explicit_limit = sim.dx / c
    dt = 4.0 * dt_explicit_limit                # 4x the explicit limit
    for _ in range(int(0.6 / dt)):
        sim.step(dt)
    print(f"max|η| = {np.max(np.abs(sim.eta)):.4f}  (still bounded)")

An explicit step diverges at $\Delta t = 4 \Delta t_\text{ex}$ . Casulli finishes the interval of length 0.6 on that same step. Phase lags. Amplitude dies faster as θ approaches 1. θ=0.5 splits the weights in half and is second-order. Past 0.6 the scheme is first-order. The code default, 0.6, already stands on that first-order side.

Crank up the time step#

Length 1 and 200 cells give $\Delta x = 0.005$ . With $c=\sqrt{9.81}\approx 3.13$ m/s the explicit limit is about 0.0016 s, and four times that is about 0.0064 s. Drop the fraction of 0.6 divided by that step and 93 steps remain. The printout only says max|η| stayed finite. It does not say how far the initial bump of 0.05 decayed.

That 0.05 is 5% of the depth $H=1$ . The linearization $\partial_t \eta + H\partial_x u = 0$ closes only while $\eta$ stays well below $H$ . Celerity is then one number, $\sqrt{gH}$ , and the tridiagonal $k$ is the same in every cell. Raise the bump to the depth and $c$ varies from place to place. This constant $k$ no longer closes the wave equation.

At θ=0.5 the same fourfold step gives $k=16\times 0.25=4$ . The diagonal surplus is still 1. The stability edge is not the matrix. It is the time weight $(θ-1/2)$ . At zero, numerical viscosity switches off and grid modes shuttle without decay.

What's missing#

Not a universal solver.

Time order changes with the knob. θ=0.5 is second-order, but it sits on the stability edge and gravity-wave fronts bend. Past 0.6 the order is one. If the gauge amplitude is the product, shrink the step or move to a higher-order IMEX scheme.

Advection stays explicit. Plain upwind is sticky. That is why an Eulerian–Lagrangian trace walks the trajectory backward and interpolates. The step can cross more than one cell and still hold, but the result moves with interpolation order and with the error in the departure point.

Hydrostatic balance $\partial p/\partial z = -\rho g$ is assumed. Breakers and drops, anywhere vertical acceleration is large, need a non-hydrostatic correction.

One obtuse triangle and the five-minute step is off the table#

Five minutes at θ=0.6 buys a 1,500-fold cut in step count and pays first-order damping of the wave height. If the tide record's amplitude is the product, do not push θ toward 1. The diagonal surplus is already 1, so a more implicit step does not make the matrix safer. The wave is what dies.

The difference equals $\partial \eta / \partial n$ only while the line of centers is perpendicular to the face. One obtuse triangle pushes the circumcenter outside the cell and that identity ends. A tangential piece remains. Whether that piece can be put back the way a diffusion flux is put back is the question asked by the diffusion fix on a non-orthogonal face.