💻 computer science

High-order adaptive discontinuous finite elements for the shallow water equations with sub-grid irregular bathymetry

This paper presents a well-balanced, mass-conserving, and positivity-preserving high-order adaptive Discontinuous Galerkin method implemented in the deal.II library that effectively simulates shallow water flows over irregular, sub-grid bathymetry and handles wet-dry fronts without requiring regularity assumptions.

Original authors: Luca Arpaia, Giuseppe Orlando, Christian Ferrarin, Luca Bonaventura

Published 2026-05-21

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Luca Arpaia, Giuseppe Orlando, Christian Ferrarin, Luca Bonaventura

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to predict how water moves through a coastal lagoon, like the famous Venice Lagoon. To do this, you need a computer model that solves the "shallow water equations"—basically, a set of rules describing how water flows, rises, and falls.

The challenge is that the ocean floor (the bathymetry) is incredibly complex. It has jagged rocks, sudden drop-offs, and hidden channels. In the past, computer models had to smooth these details out or use very simple shapes to represent the floor, which often led to inaccurate predictions, especially when water levels dropped low enough to expose dry land (wetting and drying).

This paper introduces a new, high-tech "digital microscope" for water modeling. Here is how it works, broken down into simple concepts:

1. The "High-Resolution Map" Problem

Think of the computer's grid as a giant checkerboard laid over the ocean.

The Old Way: If the ocean floor had a sharp cliff, but your checkerboard squares were huge, the computer would have to guess what the cliff looked like inside that square. It would either smooth it out (losing detail) or try to fit a wobbly curve that might wiggle wildly and cause the math to crash.
The New Way: This method allows the computer to look at the actual high-resolution map of the ocean floor, even if the map is jagged and irregular. It doesn't try to force the floor into a smooth shape. Instead, it samples the real depth at many tiny points inside each checkerboard square. It's like having a high-definition photo of the floor that the computer can "see" even if the grid squares are large.

2. The "Magic Trick" of Wet and Dry Land

One of the hardest things to simulate is when the tide goes out, leaving some areas dry and others wet.

The Problem: If the water level drops below the floor, the math gets messy. The computer might try to calculate the speed of water in a place where there is no water, which leads to "division by zero" errors or negative water depths (which is physically impossible).
The Solution: The authors created a smart switch. When the water gets too shallow in a specific area, the computer automatically lowers its "intelligence" for that spot. It switches from a complex, high-order calculation to a simple, robust "first-order" calculation (like a safety net). This ensures the water depth never goes negative and the simulation doesn't crash, even as the tide recedes and exposes the mudflats.

3. Keeping the Balance (The "Lake at Rest" Test)

A good model must be able to sit still. If you have a calm lake with a flat surface and a bumpy bottom, the water shouldn't start sloshing around just because the computer is looking at it.

The Claim: This new method is "well-balanced." It proves mathematically that if the water is supposed to be still, the computer will keep it still, even if the ocean floor underneath is jagged and irregular. It doesn't create fake waves out of thin air.

4. The "Adaptive Mesh" (The Zoom Lens)

Running a simulation over a whole lagoon is expensive for a computer. You don't need to look at every single pebble in the deep ocean, but you do need to look closely at the narrow channels where the water rushes fast.

The Analogy: Imagine a camera with a zoom lens. The new method uses Adaptive Mesh Refinement (AMR). It keeps the grid coarse (zoomed out) in calm, open areas to save time. But when it detects a fast-moving current, a sharp wave, or a tracer (like a dye) spreading out, it instantly "zooms in" (refines the grid) to capture the details.
The Result: This makes the simulation much faster. The paper shows that by using this zooming technique, they could get results just as accurate as a super-detailed, static grid, but using far less computer power and time.

5. Tracking the "Dye" (Tracers)

Often, scientists want to track things moving with the water, like pollution or salt.

The Issue: If the math for the water flow and the math for the dye aren't perfectly synchronized, the dye might magically appear or disappear, or create weird spikes.
The Fix: The authors designed the math so that the way they track the dye is perfectly consistent with how they track the water. If the water moves a certain amount, the dye moves exactly the same amount, preserving the total mass.

Summary of the Achievement

The paper presents a robust, high-order mathematical method that:

Uses real, jagged ocean floor data without smoothing it out.
Handles wet and dry land without crashing or creating fake water.
Saves time by automatically zooming in only where it's needed.
Is mathematically proven to be stable, accurate, and mass-conserving.

They tested this on idealized math problems and a realistic simulation of the Venice Lagoon. The results showed that their method could handle the complex, messy reality of coastal waters better and faster than previous approaches, capturing small-scale details that other models would miss.

Technical Summary: High-order adaptive discontinuous finite elements for the shallow water equations with sub-grid irregular bathymetry

Problem Statement
The numerical modeling of coastal flows using the shallow water equations (SWE) faces a significant challenge: the disparity between the high resolution of available bathymetric data and the coarser meshes typically employed in high-order finite element methods. Standard high-order approaches often require the bathymetry to be approximated by polynomials of the same order as the solution variables. For irregular bathymetries, this necessitates smoothing or Total Variation Diminishing (TVD) limiting to avoid Gibbs phenomena, which discards valuable high-resolution data. Furthermore, handling wetting and drying (flooding) on such meshes while maintaining mass conservation, positivity of water depth, and the "lake-at-rest" steady state (C-property) is mathematically complex, particularly when the mesh is not aligned with bathymetric gradients or wet-dry fronts.

Methodology
The authors propose a high-order Discontinuous Galerkin (DG) method implemented within the deal.II library, designed to utilize realistic, irregular bathymetry data without regularity assumptions. The core methodology involves:

Prognostic Variables: The free-surface elevation ( $\zeta$ ) is chosen as the primary prognostic variable for the continuity equation, rather than water depth ( $h$ ). This facilitates a smooth representation suitable for high-order finite elements. Water depth is recovered non-linearly at quadrature nodes via $h = \max(\zeta + z_b, 0)$ , where $z_b$ is the bathymetry sampled directly from reference data.
Sub-grid Bathymetry Representation: The bathymetry is evaluated directly at quadrature nodes without polynomial approximation or post-processing. This allows the method to recover accurate sub-grid bathymetric features even on coarse meshes.
Discretization Scheme:
- Spatial: A nodal DG method is used. The non-conservative form of the hydrostatic pressure gradient is discretized using a specific formulation (inspired by [67]) combined with hydrostatic reconstruction at element interfaces. This ensures the scheme is well-balanced for any bathymetry function.
- Temporal: An Implicit-Explicit (IMEX) Runge-Kutta (RK) scheme is employed. Friction terms are treated implicitly (using TR-BDF2) to handle stiffness in shallow regions, while advection and pressure gradients are treated explicitly (using a three-stage second-order RK scheme).
- Wetting and Drying: A degree adaptation strategy is employed. In elements where the water depth is zero or near zero (wet-dry fronts), the polynomial order is locally reduced to zero ( $r=0$ ), effectively switching to a first-order cell-centered Finite Volume approximation. This ensures stability and positivity without requiring ad hoc limiting of the bathymetry.
- Tracer Transport: A consistent conservative discretization is proposed for passive tracers. By evolving the depth-averaged concentration and ensuring the discrete forms of the continuity and tracer equations are consistent, the method avoids spurious monotonicity violations often seen in irregular bathymetry scenarios.
Adaptive Mesh Refinement (AMR): The method utilizes deal.II's capabilities for dynamic $h$ - and $p$ -adaptivity. Mesh refinement is driven by error indicators (e.g., vorticity or tracer gradients), while polynomial degree adaptation handles wet/dry regions.

Key Contributions

High-Order Sub-grid Method: The paper presents the first high-order method that simultaneously combines well-balancing, mass conservation, and positivity preservation for the SWE without assuming bathymetry regularity or requiring smoothing/limiting of the bathymetric data.
Mathematical Properties: The authors prove that the proposed scheme is:
- Well-balanced: It exactly preserves the lake-at-rest state for arbitrary bathymetries, including discontinuous ones.
- Positivity-preserving: Under a mild CFL condition dependent on the quadrature resolution, the water depth remains non-negative.
- Consistent: The discretization of tracer transport is consistent with the continuity equation (CWC property), preventing monotonicity violations.
Robustness to Irregularity: The method demonstrates robustness in the presence of irregular bathymetries and under-resolved features at the grid scale, avoiding the oscillations and instabilities common in polynomial-based bathymetry approximations.

Results
The method was validated through a series of idealized and realistic benchmarks:

Lake-at-Rest: The scheme preserved hydrostatic equilibrium up to machine precision on both Cartesian and distorted meshes with discontinuous bathymetries, even with wet-dry interfaces.
Channel Flow: In tests with unresolved bathymetric obstacles, the proposed sub-grid approach yielded accurate steady-state solutions. In contrast, polynomial approximations of the bathymetry exhibited spurious oscillations or failed to converge to the correct steady state.
Moving Wet-Dry Fronts (Thacker Test): The method correctly handled periodic oscillations with wet-dry transitions, maintaining mass conservation and tracer consistency within numerical round-off errors.
Travelling Vortex: Adaptive Mesh Refinement (AMR) simulations achieved the same accuracy as static high-resolution meshes but with significantly fewer degrees of freedom (up to 9x reduction) and reduced computational time (5–7x faster).
Venice Lagoon Benchmark: A realistic simulation of tidal circulation in the Venice Lagoon demonstrated the method's ability to resolve complex geometries and sub-grid channels. The sub-grid approach produced stronger velocities and more accurate flooding patterns compared to limited piecewise linear bathymetry approximations. While dynamic AMR introduced a small mass conservation error ( $\sim 10^{-5}$ ), this was deemed acceptable relative to typical boundary data uncertainties.

Significance and Claims
The paper claims that this work provides a robust and accurate framework for coastal modeling that fully exploits high-resolution bathymetric datasets. By decoupling the bathymetry representation from the polynomial basis of the solution, the method avoids the trade-off between accuracy and stability inherent in traditional high-order schemes. The authors emphasize that their approach allows for the simulation of complex coastal environments with irregular bathymetries and wetting/drying fronts without the need for data regularization. The integration with deal.II enables efficient parallelization and dynamic adaptivity, making the method suitable for high-fidelity, computationally efficient coastal flow simulations. The work serves as a foundation for future developments, such as implicit treatment of pressure gradients to relax time-step restrictions based on celerity.