Mathematical Background

This page outlines the mathematics underpinning ANUGA. It is a web-searchable version of Chapter 16 of the ANUGA User Manual. ANUGA has been applied to hydrodynamic modelling of coastal inundation [Nielsen2005] and many other shallow water problems.

Shallow Water Wave Equations

The shallow water wave equations are a system of differential conservation equations describing the flow of a thin layer of fluid over terrain:

\[\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{E}}{\partial x} + \frac{\partial \mathbf{G}}{\partial y} = \mathbf{S}\]

The vector of conserved quantities is

\[\begin{split}\mathbf{U} = \begin{bmatrix} h \\ uh \\ vh \end{bmatrix}\end{split}\]

where \(h\) is water depth, \(uh\) is momentum in the \(x\)-direction, and \(vh\) is momentum in the \(y\)-direction. Other quantities entering the system are bed elevation \(z\) and stage (absolute water level) \(w\), related by \(w = z + h\) at all times.

The flux vectors in the \(x\) and \(y\) directions are

\[\begin{split}\mathbf{E} = \begin{bmatrix} uh \\ u^2 h + \tfrac{1}{2}gh^2 \\ uvh \end{bmatrix}, \qquad \mathbf{G} = \begin{bmatrix} vh \\ uvh \\ v^2 h + \tfrac{1}{2}gh^2 \end{bmatrix}\end{split}\]

and the source term (gravity and friction) is

\[\begin{split}\mathbf{S} = \begin{bmatrix} 0 \\ -gh(z_x + S_{fx}) \\ -gh(z_y + S_{fy}) \end{bmatrix}\end{split}\]

where \(S_f\) is the bed friction modelled by Manning’s resistance law:

\[S_{fx} = \frac{u \eta^2 \sqrt{u^2 + v^2}}{h^{4/3}}, \qquad S_{fy} = \frac{v \eta^2 \sqrt{u^2 + v^2}}{h^{4/3}}\]

with \(\eta\) the Manning resistance coefficient. The model does not include kinematic viscosity or dispersion (though Kinematic_viscosity_operator can be added as an explicit operator).

Finite Volume Method

ANUGA uses a finite-volume method to solve the shallow water wave equations [ZR1999]. The study area is represented by a mesh of triangular cells; the conserved quantities \(h\), \(uh\), \(vh\) are stored at the centroid of each cell. The size of the triangles may vary to provide greater resolution in regions of particular interest.

Cell update equation

Integrating the conservation equations over triangular cell \(T_i\) and applying the divergence theorem gives the rate equation for the cell-averaged conserved quantities:

\[\frac{d\mathbf{U}_i}{dt} + \frac{1}{A_i} \sum_{j \in \mathcal{N}_i} \mathbf{H}_{ij}\, l_{ij} = \mathbf{S}_i\]

where

• \(\mathbf{U}_i\) — conserved quantities averaged over cell \(i\)

• \(A_i\) — area of cell \(i\)

• \(\mathcal{N}_i\) — set of indices \(j\) of cells neighbouring cell \(i\)

• \(l_{ij}\) — length of edge \(e_{ij}\) between cells \(i\) and \(j\)

• \(\mathbf{H}_{ij}\,l_{ij}\) — outward normal flux across edge \(e_{ij}\)

• \(\mathbf{S}_i\) — average source term over cell \(i\)

Reconstruction

A second-order piecewise-linear reconstruction of the conserved quantities is built within each cell from the cell-average values of that cell and its neighbours. The slope is limited (see Slope Limiting below) to suppress spurious oscillations. The reconstruction is allowed to be discontinuous across edges. The values on either side of edge \(e_{ij}\) are denoted

\[\mathbf{U}^i_{ij} = \lim_{\mathbf{x} \to m_{ij}} \mathbf{U}_i(\mathbf{x}), \qquad \mathbf{U}^j_{ij} = \lim_{\mathbf{x} \to m_{ij}} \mathbf{U}_j(\mathbf{x})\]

where \(m_{ij}\) is the midpoint of edge \(e_{ij}\).

Numerical flux — central-upwind scheme

The normal flux \(\mathbf{H}(\mathbf{U}) = \mathbf{E}(\mathbf{U})\,n_1 + \mathbf{G}(\mathbf{U})\,n_2\) is approximated using the central-upwind scheme of Kurganov, Noelle & Petrova [KurNP2001]:

\[\mathbf{H}_{ij} = \frac{a^+_{ij}\,\mathbf{H}(\mathbf{U}^i_{ij}) - a^-_{ij}\,\mathbf{H}(\mathbf{U}^j_{ij})} {a^+_{ij} - a^-_{ij}} + \frac{a^+_{ij}\,a^-_{ij}}{a^+_{ij} - a^-_{ij}} \bigl[\mathbf{U}^j_{ij} - \mathbf{U}^i_{ij}\bigr]\]

with one-sided wave speed bounds

\[\begin{split}a^+_{ij} &= \max\!\left\{ u^i_{ij} + \sqrt{g h^i_{ij}},\; u^j_{ij} + \sqrt{g h^j_{ij}},\; 0 \right\} \\[4pt] a^-_{ij} &= \min\!\left\{ u^i_{ij} - \sqrt{g h^i_{ij}},\; u^j_{ij} - \sqrt{g h^j_{ij}},\; 0 \right\}\end{split}\]

Time stepping

An explicit Euler scheme with adaptive timestep is used:

\[\mathbf{U}^{n+1}_i = \mathbf{U}^n_i - \frac{\Delta t}{A_i} \sum_{j \in \mathcal{N}_i} \mathbf{H}_{ij}\,l_{ij} + \Delta t\,\mathbf{S}_i\]

Stability requires the timestep to satisfy the Courant–Friedrichs–Lewy (CFL) condition:

\[\Delta t \le \min_{i}\; \min_{j \in \mathcal{N}_i} \left( \frac{r_i}{a_{ij}},\; \frac{r_j}{a_{ij}} \right)\]

where \(a_{ij} = \max\{a^+_{ij},\, -a^-_{ij}\}\) and \(r_i\) is the inradius (inscribed circle radius) of cell \(i\).

Flux Limiting

Near wet/dry boundaries the water depth \(h\) becomes very small, making the velocity recovery \(u = uh/h\) numerically unreliable and potentially producing unphysical speeds.

ANUGA replaces the raw velocity calculation with the limited approximation

\[\hat{u} = \frac{uh \cdot h}{h^2 + h_0}, \qquad \hat{v} = \frac{vh \cdot h}{h^2 + h_0}\]

where \(h_0\) is a regularisation parameter. Taking limits:

\[\lim_{h \to 0} \hat{u} = 0, \qquad \lim_{h \to \infty} \hat{u} = \frac{uh}{h} = u\]

so the limiter smoothly recovers zero velocity at dryness and the true velocity in deep water.

ANUGA exposes a global minimum-depth parameter \(H_0\) (typically \(10^{-3}\) m, set in anuga/config.py as minimum_allowed_height). Setting

\[h_0 = H_0^2\]

provides a good balance between accuracy and stability. At depth \(h = N H_0\) the predicted speed is scaled by

\[\frac{1}{1 + 1/N^2}\]

which converges quadratically to the true value with \(N\). The limiter is applied only for depths less than \(10\,H_0\) (roughly 1 cm for the default \(H_0\)), where it affects the computed velocity by less than 1 %.

Slope Limiting

A multidimensional slope-limiting technique (MinMod limiter) achieves second-order spatial accuracy and prevents spurious oscillations. Near the bed, the limiter must additionally ensure that no negative depths occur.

Let \(w\), \(z\), \(h\) be stage, bed elevation, and depth at the centroid, and let \(w_i\), \(z_i\), \(h_i\) be the corresponding vertex values. Define the minimum vertex depth as \(h_{\min} = \min_i h_i\).

Two reconstructed stages are computed:

• \(\tilde{w}_i\) — stage from a gradient limiter applied to stage only (suitable for deep water where the bed slope can be ignored)

• \(\bar{w}_i\) — stage from a gradient limiter applied to depth, respecting the bed slope (suitable for shallow water near wet/dry fronts)

The balanced stage is the linear combination

\[w_i = \alpha\,\tilde{w}_i + (1-\alpha)\,\bar{w}_i, \qquad \alpha \in [0,1]\]

In deep water (\(h_{\min} \ge \varepsilon\)) we set \(\alpha = 1\) (stage-limited reconstruction only).

In shallow water (\(h_{\min} < \varepsilon\)) we choose the largest \(\alpha\) that still guarantees non-negative vertex depths:

\[\alpha = \min_i \frac{\bar{h}_i - \varepsilon}{\bar{h}_i - \tilde{h}_i}\]

where \(\bar{h}_i = \bar{w}_i - z_i\) and \(\tilde{h}_i = \tilde{w}_i - z_i\). Should \(\alpha\) fall outside \([0, 1]\) it is clamped to that interval.

See also

ANUGA User Manual — Chapter 16: Mathematical Background is the original source for this page and includes the parallel implementation theory, scalability results, and the full TikZ mesh diagrams.

References

[ZR1999]

Zoppou, C. & Roberts, S. (1999). Catastrophic collapse of water supply reservoirs in urban areas. Journal of Hydraulic Engineering, 125(7), 686–695.

[KurNP2001]

Kurganov, A., Noelle, S. & Petrova, G. (2001). Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton–Jacobi equations. SIAM Journal on Scientific Computing, 23(3), 707–740.

[Nielsen2005]

Nielsen, O., Roberts, S., Gray, D., McPherson, A. & Hitchman, A. (2005). Hydrodynamic modelling of coastal inundation. MODSIM 2005 International Congress on Modelling and Simulation, 518–523.