Source-linked AI summary
Thetis coastal ocean model: discontinuous Galerkin discretization for the three-dimensional hydrostatic equations
Tuomas Kärnä, Stephan C. Kramer, Lawrence Mitchell, David A. Ham, Matthew D. Piggott, António M. Baptista
TL;DR
Coastal unstructured-grid models can be too diffusive or computationally expensive, limiting their use for multiscale coastal flows. The paper presents a DG discretization for the hydrostatic equations and evaluates its conservation and accuracy. The tests report conservation of volume and tracers and second-order accuracy under barotropic forcing.
Problem
Unstructured-grid coastal models can suffer excessive numerical diffusion or high computational cost, while numerical mixing can smear fronts and misrepresent coastal mixing.
Source: §1 · p. 2Method
The paper presents a discontinuous Galerkin discretization of the three-dimensional hydrostatic equations in Thetis, an open-source coastal ocean circulation model.
Source: §1 · pp. 3–4Results
The standing-wave test reports volume and tracer conservation, negligible temperature overshoots with slope limiters, and second-order accuracy under barotropic forcing.
Source: §5.1 · pp. 16–17Takeaways & Limitations
The standing-wave results verify that the solver behaves as expected and yields second-order accuracy under barotropic forcing.
Source: §5.1 · p. 17The chosen function spaces are not mimetic and do not exactly conserve all continuous-equation properties, such as enstrophy; a mimetic alternative increased computational cost without corresponding accuracy improvement in preliminary experiments.
Source: §3.1 · p. 8
Abstract
from arXiv · showhide
Unstructured grid ocean models are advantageous for simulating the coastal ocean and river-estuary-plume systems. However, unstructured grid models tend to be diffusive and/or computationally expensive which limits their applicability to real life problems. In this paper, we describe a novel discontinuous Galerkin (DG) finite element discretization for the hydrostatic equations. The formulation is fully conservative and second-order accurate in space and time. Monotonicity of the advection scheme is ensured by using a strong stability preserving time integration method and slope limiters. Compared to previous DG models advantages include a more accurate mode splitting method, revised viscosity formulation, and new second-order time integration scheme. We demonstrate that the model is capable of simulating baroclinic flows in the eddying regime with a suite of test cases. Numerical dissipation is well-controlled, being comparable or lower than in existing state-of-the-art structured grid models.
1 Introduction
Coastal ocean models must resolve multiscale dynamics, but unstructured-grid models can be overly diffusive or computationally expensive. The paper proposes a DG hydrostatic solver with design criteria and implementation choices aimed at accuracy, conservation, and efficient simulation.
1 Introduction
- The proposed DG hydrostatic solver targets unstructured grids, where excessive numerical diffusion and high computational cost are recognized drawbacks.Unstructured grids can represent complex coastal topography and local features, while structured grids offer better computational performance.Source: §1 · pp. 2–3
- Numerical mixing can smear coastal fronts and misrepresent saline intrusion or river-plume mixing, motivating attention to accuracy in these strongly structured flows.Coastal domains include strong density gradients, fronts, internal waves, hydraulic jumps, and baroclinic eddies.Source: §1 · p. 2
- The design criteria combine a vertically extruded layered mesh, accurate free-surface dynamics, second-order monotone tracer advection, explicit 3D time integration except vertical diffusion, and low numerical mixing.The tracer scheme uses linear DG elements, slope limiting, and SSP time integration; an ALE formulation accounts for free-surface motion and supports generic vertical grids.Source: §1 · p. 3
- The paper focuses on numerical closures, using viscosity schemes mainly for stability rather than to represent unresolved physical mixing.The discussion contrasts Laplacian viscosity, which can dissipate larger scales, with biharmonic viscosity, which more strongly dissipates smaller scales.Source: §1 · pp. 3–4
- Thetis implements both depth-averaged 2D circulation and full 3D hydrostatic models, with the paper focusing on the 3D model.The implementation uses Firedrake for extruded-mesh support, flexible weak-form descriptions, and automated code generation.Source: §1 · p. 4
2 Governing equations
The governing system couples hydrostatic horizontal momentum, tracer transport, diagnosed vertical velocity, and free-surface dynamics. It separates depth-averaged and deviation modes, then specifies their coupling, density calculation, and viscosity closures.
2 Governing equations
- Horizontal momentum combines Coriolis, pressure-gradient, and viscous terms, while temperature and salinity follow advection-diffusion equations and vertical velocity is diagnosed from continuity.The model also applies quadratic bottom stress, treated implicitly with the bottom-element velocity magnitude fixed at its old value.Source: §2 · p. 5
2.1 Mode splitting
- The depth-averaged system carries fast rotational surface gravity waves; advection and viscosity remain unsplit in the deviation equation under a long-time-step assumption.The Coriolis term contains only slow modes, and vertical velocity needs no separate split because it appears only in the unsplit advection term.Source: §2.1 · p. 6
2.2 Coupling 2D and 3D modes
- A vertically constant forcing G corrects the deviation velocity's nonzero depth average, so the 2D and 3D equations sum to the full momentum equation.The correction enforces zero depth average for the deviation velocity.Source: §2.2 · p. 6
2.3 Equation of state
- The test cases use a linear equation of state with salinity excluded from density, while Thetis also implements a nonlinear equation of state.The linear form allows density to be computed locally at tracer-field nodes.Source: §2.3 · p. 7
2.4 Viscosity and turbulence closure
- The model uses constant Laplacian horizontal viscosity; vertical viscosity is usually constant but sometimes follows a gradient-Richardson-number parametrization.More sophisticated turbulence closures are deferred to future work.Source: §2.4 · p. 7
3 Finite element discretization
This section assigns finite-element spaces and derives weak forms for the coupled system. Its design addresses tracer monotonicity and momentum stability, while acknowledging that the chosen spaces are not mimetic.
3.1 Function spaces
- Tracer consistency is supported by placing vertical velocity in the tracer space, but the chosen spaces do not conserve enstrophy exactly and mimetic alternatives may raise computational cost without improving accuracy.Tracer consistency means a constant tracer remains free of spurious local extrema; the discrete tracer equation must reduce to continuity for a constant tracer.Source: §3.1 · p. 8
3.2 2D system
- The 2D weak form uses approximate Riemann-solver interface values, with a linear Roe solution specified for the free-surface variable.The divergence and external pressure-gradient terms are integrated by parts.Source: §3.2 · p. 9
3.3 Momentum equation
- The 3D momentum weak form integrates advection and viscosity by parts and adds a Lax–Friedrichs flux to the internal pressure gradient to reduce velocity noise and spurious mixing.The flux parameter is given as γlf = |u|.Source: §3.3 · p. 10
3.4 Tracer equation
The discretization treats tracer transport, diffusion, continuity, and internal pressure through weak formulations and stabilization. Slope limiters are applied after advection updates to support positivity.
3.4 Tracer equation
- The tracer weak formulation is derived analogously to the momentum formulation, but it does not use the Lax–Friedrichs flux.The flux distinction is specific to the tracer equation.Source: §3.4 · p. 10
3.5 Symmetric Interior Penalty stabilization
- The discretization is unstable for elliptic operators, so tracer diffusion uses Symmetric Interior Penalty Galerkin stabilization, with the viscosity terms formulated similarly.The penalty factor depends on mesh quality and local element length scale; the paper uses γ = 5.Source: §3.5 · pp. 10–11
3.6 Continuity equation
- Vertical velocity is computed diagnostically from continuity using a weak form; bottom-surface terms vanish under the impermeability constraint.Source: §3.6 · p. 11
3.7 Computing the internal pressure gradient
- Density is computed diagnostically from the equation of state, and the baroclinic head is obtained by vertically integrating density anomaly and solving weakly with boundary conditions.The paper uses a linear equation of state in this work; the free-surface terms vanish because the baroclinic head is zero there by definition.Source: §3.7 · pp. 11–12
3.8 Slope limiters
- Three-dimensional slope limiters are applied to tracer and horizontal-velocity fields after each advection-operator update to ensure positivity.Source: §3.8 · p. 12
4 Time integration
The coupled 2D–3D solver uses ALE mesh updates to account for free-surface movement, with a second-order time integration scheme. Its time-step constraints are set by explicit stability limits.
4.1 ALE mesh formulation
- ALE mesh movement accounts for free-surface changes by uniformly stretching or compressing vertical coordinates while maintaining a conforming mesh.The discontinuous elevation field is projected to a continuous space; the resulting discrepancy is typically negligible because elevation jumps are small.Source: §4.1 · p. 12
4.2 Coupled time integration scheme
- The 2D system is solved first to determine elevation and mesh geometry, enabling mesh-velocity calculation before the corresponding 3D ALE updates.The explicit 3D equations use two SSP stages, with slope limiting and coupling corrections during the coupled update.Source: §4.2 · pp. 13–15
- The coupled scheme is second-order for all terms, combining SSPRK(2,2) updates with implicit treatment of selected vertical and bottom-friction terms.Only vertical diffusion, vertical viscosity, and bottom friction are implicit; the 2D scheme is implicit and compatible with SSPRK(2,2).Source: §4.2 · pp. 13–15
4.3 Choosing the time step
- The coupled system’s stability follows the explicit SSP scheme’s conditions; horizontal advection, vertical advection, and explicit viscosity impose additional time-step constraints.The simulations use a constant time step selected a priori from mesh-based velocity and viscosity scales.Source: §4.3 · pp. 15–16
5 Test cases
The test suite checks conservation and convergence under barotropic forcing, then verifies discrete baroclinic operators using manufactured solutions. The tests report conserved volume and tracers and expected convergence behavior.
5 Test cases
- The test suite progresses from a barotropic standing wave to baroclinic manufactured solutions, lock exchange, a rotating eddy, and DOME overflow.The standing-wave setup uses a closed rectangular channel with an initial sinusoidal elevation perturbation and passive salinity and temperature tracers.Source: §5 · p. 16, §5.1 · p. 16
5.1 Standing wave
- Volume was conserved to O(10^-15), 2D volume to O(10^-16), and tracer mass to O(10^-12); limiters reduced temperature overshoots from O(10^-2) to negligible levels.Constant salinity deviated by only O(10^-9), and tracer consistency supported the integrity of the ALE formulation.Source: §5.1 · p. 16
- The standing-wave test showed second-order velocity convergence and elevation convergence at a rate of 3.2, consistent with expected solver accuracy under barotropic forcing.Figure 1 presents L2 error convergence across element sizes from 3000 to 300.Source: §5.1 · p. 17, Figure 1 · p. 17
Appendix A.
The tests report second-order convergence, controlled mixing in lock exchange and eddying flows, and a DOME plume whose structure agrees qualitatively with prior studies. Results also characterize solver scaling and note that Thetis performance is not fully optimized.
Appendix A.
- Across the MMS refinements, prognostic variables show second-order convergence, while diagnostic vertical velocity converges linearly; tested viscosity, diffusivity, and bottom-friction operators also show second-order convergence.The reported convergence supports the discretization of advection, pressure-gradient, and Coriolis terms.Source: Appendix A · p. 18, Figure 2 · p. 19
5.3 Lock exchange
- Lock-exchange front propagation agrees with theory, and at Reh = 250 temperature overshoots are about 10−5 °C with slope limiters, versus up to 30 °C without them.The cited results attribute near-monotone tracer advection to SSP time integration and slope limiters.Source: §5.3 · pp. 19–20, Figure 4 · p. 21
- The lock-exchange RPE results indicate mixing-related energy loss of roughly 5%, similar to reported results for ROMS, MITgcm, and MOM.At 17 h, normalized RPE values across four simulations are 0.612, 1.13, 2.35, and 3.11 ×10−5.Source: §5.3 · pp. 20–21
5.4 Baroclinic eddies
- In the eddying channel, lowering horizontal viscosity yields more mesoscale features, while normalized RPE is roughly twice as high at Reh = 20 as at Reh = 2.RPE-change rates also increase with grid Reynolds number and with coarser mesh resolution.Source: Figure 5 · p. 22, §5.4 · p. 23, Figure 6 · p. 24
- Thetis's RPE-change rates at Reh = 20 are 4.3 × 10−4 and 2.2 × 10−4 W m−2 for 10 and 4 km meshes, below reported MITgcm, MOM, and POP values.On a hexagonal mesh, MPAS-Ocean reports smaller rates; on a quad mesh its values are close to Thetis performance.Source: §5.4 · pp. 23–24, Figure 6 · p. 24
- Thetis runs 3.3× faster than SLIM 3D on one core, 4.0× faster on 24 cores, and 2.2× faster on 144 cores in the reported comparison.The models use identical setup, mesh, and time step; Thetis performance is not fully optimized.Source: §5.4 · p. 25
5.5 DOME
- The DOME simulation forms a roughly 150 km coastal plume with eddies and internal waves; most tracer remains within 200 m of the bottom.The plume follows the bottom slope, with its most vigorous eddies in the first 300 km after the inlet.Source: §5.5 · p. 26, Figure 9 · p. 27
- The DOME plume's density changes as it travels downstream: inflowing water is entrained into higher density-spectrum levels, and the mean density anomaly ranges from 0.8 to 1.5 kg m−3.The anomaly peaks around 1.28 kg m−3; the reported rate of entrainment is used as a mixing metric.Source: §5.5 · pp. 27–28
6 Conclusions
The DG coastal-ocean solver is second-order accurate, fully conservative, and preserves monotonicity. Test cases reproduce expected behavior, including baroclinic eddies, with numerical mixing comparable to established models; several extensions remain future work.
6 Conclusions
- Test cases reproduce expected physical behavior, including baroclinic eddies, with numerical mixing comparable to MITgcm, ROMS, and MPAS-Ocean.The paper describes this as a contribution toward more accurate next-generation coastal ocean models.Source: §6 · p. 28
- The DG solver is second-order accurate in space and time, fully conservative, and preserves monotonicity.Source: §6 · p. 28
- Future work includes spherical grids, biharmonic viscosity, turbulence closures, wetting-drying, an adjoint solver, and improved computational efficiency and parallel scaling.Source: §6 · p. 29
7 Code availability
The experiment code is publicly available, and the paper identifies archived versions of Firedrake components and Thetis for reproducibility.
7 Code availability
- All code used for the experiments is publicly available.Source: §7 · p. 29
- Firedrake and its components are available from firedrakeproject.org, and Thetis is available from thetisproject.org.Source: §7 · p. 29
- Exact software versions used to produce the results are archived, including Firedrake components and Thetis.The paper points to Zenodo records for both.Source: §7 · p. 29
8 Data availability
The manuscript reports that it used no external data.
8 Data availability
- No external data were used in the manuscript.Source: §8 · p. 29
Appendix A: Source terms for the baroclinic MMS test
For the baroclinic MMS test, analytical velocity and temperature fields provide the basis for deriving steady-state fields and evaluating equation terms. Those terms enter the equations as source terms in the weak form.
Appendix A: Source terms for the baroclinic MMS test
- Analytical velocity and temperature fields are used to derive the steady-state solution for the remaining fields.Source: Appendix A · p. 29
- The different terms in the momentum and tracer equations are evaluated for the test.Source: Appendix A · p. 29
- The evaluated terms are added as right-hand-side source terms and, in weak form, multiplied by the test function and integrated over the domain.The solutions were derived using SymPy.Source: Appendix A · p. 30
Appendix B: CPU cost comparison against SLIM
In the baroclinic eddies test case, Thetis ran faster than SLIM across the tested core counts. The comparison used identical mesh and simulation settings.
Appendix B: CPU cost comparison against SLIM
- Thetis ran approximately 3.3× faster than SLIM on one core, 4.0× faster on 24 cores, and 2.2× faster on 144 cores.Table B1 reports CPU time for the baroclinic eddies test case.Source: Appendix B · p. 31, Table B1 · p. 31
- The scaling comparison defines speed-up as T0/Ti and parallel efficiency as speed-up divided by core count; ideal efficiency is unity.Tests ran on a Linux cluster with 16-core Intel Xeon E5620 processors and Mellanox Infiniband interconnect.Source: Appendix B · p. 31
- Both models used an identical triangular mesh with 4 km resolution and 40 vertical levels, ν = 20 m2 s−1, and a 140 s time step.Wall-clock time was recorded over 40 iterations.Source: Table B1 · p. 31