Thetis coastal ocean model: discontinuous Galerkin discretization for the three-dimensional hydrostatic equations

Abstract

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.

Introduction

Numerical modeling of the coastal ocean is important for many environmental and industrial applications. Typical scenarios include modeling circulation at regional scales, coupled river–estuary–plume systems, river networks, lagoons, and harbors. Length scales range from some tens of meters in rivers and embayments to tens of kilometers in the coastal ocean; water depth ranges from less than a meter to kilometer scale at the shelf break. The timescales of the relevant processes range from minutes to hours, yet typical simulations span weeks or even decades. The dynamics are highly non-linear, characterized by local small-scale features such as fronts and density gradients, internal waves, and baroclinic eddies. These physical characteristics imply that coastal ocean modeling is intrinsically multi-scale, which imposes several technical challenges.

Most coastal ocean models solve the hydrostatic Navier–Stokes equations under the Boussinesq approximation – a valid approximation for mesoscale and submesoscale (1 km) processes. Small-scale processes ($<\mathrm{100}$ m) are, however, inherently three-dimensional where non-hydrostatic effects can be important, especially in areas with pronounced density structure and stratification . Non-hydrostatic modeling requires very high horizontal mesh resolution, which is currently only feasible in relatively small subregions (e.g., at the mouth of an estuary; ) due to its high computational cost.

Historically, regional ocean models have used structured, (deformed) rectilinear lattice grids. Although structured grids offer better computational performance , unstructured grids are generally preferred in coastal domains as they can better represent the complex coastal topography and local features . Due to the large geometrical aspect ratio of the oceans (length versus depth), most models utilize computational grids that are layered in the vertical direction. Typical approaches include the terrain-following sigma levels , equipotential $z$ levels , isopycnal coordinates , and their generalizations e.g.,.

In this article, we focus on solving the hydrostatic equations on an unstructured grid. While many unstructured grid models exist, their drawbacks tend to be excessive numerical diffusion that smooths out important physical features and/or high computational cost. To address these issues, we propose a novel finite element solver for the hydrostatic equations, based on discontinuous Galerkin discretization methods.

Maintaining high numerical accuracy is crucial in ocean applications. The ocean is a forced dissipative system where the mixing of water masses only takes place at the molecular level . In practice, however, the finite grid resolution and numerical schemes used by the model introduce mixing rates of tracers and momentum that can be orders of magnitude larger than physical mixing . Such spurious, numerical mixing is often dominated by the discretization of advection , but it can arise from other components as well, such as (implicit) time integration methods or various filters introduced to improve numerical stability . In addition, wetting and drying schemes may introduce additional dissipation in order to stabilize the barotropic equation in the drying regime. We reserve consideration of this important latter topic for a future publication.

In global circulation models, numerical mixing is a major bottleneck as (diapycnal) diffusion is very low in the deep ocean basins and water masses can remain largely unchanged for hundreds of years . Numerical mixing can, however, be a major issue in coastal domains as well: coastal oceans are characterized by strong density gradients, fronts between water masses (e.g., in river plumes), small-scale dynamics (e.g., internal waves and hydraulic jumps), and baroclinic eddies. An overly diffusive model can, therefore, fail to capture many essential physical features of these domains: it can smear out fronts, underestimate the intrusion of saline waters into embayments , or misrepresent mixing in river plumes.

The most common spatial discretization scheme is the finite volume (FV) method, used in the MITgcm , GETM , ROMS , MPAS-Ocean , UnTRIM , FVCOM , SUNTANS , FESOM2 , and others. The FV method is well suited for advection-dominated problems, provides strict conservation of volume and mass, and yields good computational performance. FV methods are nominally only first-order accurate, but higher-order approximations can be introduced by increasing the size of the numerical stencil (e.g., in high-order advection schemes; ).

Some unstructured grid models are based on the continuous Galerkin finite element (FE) method or hybrid FE–FV formulations. Such models include ADCIRC , SELFE , and SCHISM , and the earlier version of FESOM . The continuous FE method is ideal for solving elliptic equations but requires stabilization for advection seeand references therein. In addition, these methods involve solving a fully coupled global system which is less efficient in parallel applications compared to the FV method .

In recent years, discontinuous Galerkin (DG) methods have gained attention in geophysical modeling . DG discretization resembles the FV method because it is local (i.e., elements are only connected by inter-element fluxes), fully conservative, and well-suited for advective problems, yet it offers higher-order accuracy. This article presents a DG discretization for the hydrostatic equations. Our goal is to design an efficient unstructured grid solver where numerical accuracy is not compromised. Specifically, we aim to meet the following design criteria:

• a vertically extruded, layered mesh;

• accurate representation of free surface dynamics;

• a second-order accurate, monotone tracer advection scheme;

• explicit time integration of 3-D variables (except for vertical diffusion); and

• low numerical mixing.

Based on the advection scheme requirements, we have chosen to use linear discontinuous Galerkin elements for tracers, combined with a slope limiter and a strong stability-preserving (SSP) time integration scheme . This choice ensures that the scheme is second order in smooth areas, while slope limitation combined with the SSP time integration scheme ensure monotonicity (i.e., no overshoots). The movement of the free surface is taken into account with an arbitrary Lagrangian–Eulerian (ALE) formulation , where the mesh moves in the vertical direction. The ALE formulation guarantees strict local and global conservation of volume and tracers and allows for the use of generic vertical grids .

All numerical ocean models include some form of friction, either in the form of a numerical closure or a physical parameterization . Numerical closure involves adding a sufficient amount of dissipation to maintain numerical stability. There is a wealth of literature about stable finite volume e.g., and finite element discretizations e.g., for rotational shallow water equations. Most of these schemes are stable for external gravity waves and hence do not require any additional dissipation. Solving the 3-D hydrostatic equations under strong baroclinic forcing, however, generates noise at the grid scale that does require dampening. A common approach is to add some form of viscosity proportional to the grid Reynolds number . argue that conventional Laplacian viscosity has too wide a spectrum and tends to dissipate physically relevant (larger) scales too much. They show that biharmonic viscosity dissipates smaller scales more and is thus more appropriate for removing noise at the grid scale. In contrast to numerical closures, physical parameterizations aim to represent unresolved subgrid-scale processes, such as strong lateral mixing near coasts or mixing at the bottom boundary layers. In this article, we focus on numerical closures; the presented viscosity schemes are mostly motivated by numerical stability considerations.

In this article, we present an efficient DG implementation of the three-dimensional hydrostatic equations. The model is implemented in the Thetis project – an open-source coastal ocean circulation model freely available online (see http://thetisproject.org, last access: 25 October 2018). Thetis implements both a 2-D depth-averaged circulation model and a full 3-D hydrostatic model, the latter of which is discussed herein.

Thetis is implemented using the Firedrake finite element modeling platform https://www.firedrakeproject.org/, last access: 25 October 2018;. We have chosen Firedrake because of its flexibility and support for extruded meshes . Firedrake uses high-level abstractions for describing the weak formulation of partial differential equations, specifically the Unified Form Language , and automated code generation to produce efficient C code and just-in-time compilation. As such, it is an extremely flexible modeling framework that does not sacrifice computational efficiency; it is also an ideal platform for experimenting and benchmarking different discretizations. Automated code generation can also support different target hardware architectures, making it attractive for current and emerging high-performance computing platforms. In addition, Firedrake can automatically derive the adjoint of the forward model , permitting inverse modeling applications such as parameter optimization and data assimilation.

The governing equations are presented in Sect. , followed by their DG finite element discretization in Sect. . The second-order coupled time integration scheme is described in Sect. . Numerical tests are presented in Sect. .

Governing equations

Let $\mathrm{\Omega }$ be the three-dimensional domain that spans from the sea floor $z=-h(x$, $y)$ to the free surface $z=\mathit{\eta }(x$, $y)$; the bottom and top surfaces are denoted by ${\mathrm{\Gamma }}_{\mathrm{b}}$ and ${\mathrm{\Gamma }}_{\mathrm{s}}$, respectively. Total water column depth is thus $H=\mathit{\eta }+h$. The two-dimensional horizontal domain is denoted by ${\mathrm{\Gamma }}_{\mathrm{0}}$.

The horizontal momentum equation reads $$\begin{array}{cc}{\displaystyle \frac{\partial \mathit{u}}{\partial t}}& {\displaystyle }+{\mathbf{\nabla }}_{\mathrm{h}}\cdot \left(\mathit{u}\mathit{u}\right)+{\displaystyle \frac{\partial \left(w\mathit{u}\right)}{\partial z}}+f{\mathit{e}}_z\wedge \mathit{u}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathbf{\nabla }}_{\mathrm{h}}p\\ & {\displaystyle }& {\displaystyle }={\mathbf{\nabla }}_{\mathrm{h}}\cdot \left({\mathit{\nu }}_{\mathrm{h}}{\mathbf{\nabla }}_{\mathrm{h}}\mathit{u}\right)+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\nu }{\displaystyle \frac{\partial \mathit{u}}{\partial z}}\right),\end{array}$$ where $\mathit{u}=(u$, $v)$ and $w$ denote the horizontal and vertical velocity, respectively; ${\mathbf{\nabla }}_{\mathrm{h}}$ is the horizontal gradient operator; $\wedge$ denotes the cross product operator; $f$ is the Coriolis parameter; ${\mathit{e}}_z$ is the vertical unit vector; $p$ is the pressure; and ${\mathit{\nu }}_{\mathrm{h}}$ and $\mathit{\nu }$ are the horizontal and vertical diffusivity, respectively. Water density is defined as $\mathit{\rho }={\mathit{\rho }}_{\mathrm{0}}+{\mathit{\rho }}^{\prime }(T$, $S$, $p)$, where $T$ and $S$ stand for temperature and salinity, respectively, and ${\mathit{\rho }}_{\mathrm{0}}$ is a constant reference density.

Under the hydrostatic assumption, the horizontal pressure gradient can be written as a combination of external, internal, and atmospheric pressure gradients: $${\displaystyle }{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathbf{\nabla }}_{\mathrm{h}}p=g{\mathbf{\nabla }}_{\mathrm{h}}\mathit{\eta }+g{\mathbf{\nabla }}_{\mathrm{h}}r+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathbf{\nabla }}_{\mathrm{h}}{p}_{\mathrm{atm}},$$ where $p_{\mathrm{atm}}$ is the atmospheric pressure acting on the sea surface, and $${\displaystyle }{\displaystyle }r={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\underset{z}{\overset{\mathit{\eta }}{\int }}{\mathit{\rho }}^{\prime }\mathrm{d}{z}^{\prime }$$ is the baroclinic head. For brevity, the internal pressure gradient field is denoted as ${\mathit{F}}_{\mathrm{pg}}=g{\mathbf{\nabla }}_{\mathrm{h}}r$.

Neglecting atmospheric pressure, the full horizontal momentum equation reads $$\begin{array}{cc}{\displaystyle \frac{\partial \mathit{u}}{\partial t}}& {\displaystyle }+{\mathbf{\nabla }}_{\mathrm{h}}\cdot \left(\mathit{u}\mathit{u}\right)+{\displaystyle \frac{\partial \left(w\mathit{u}\right)}{\partial z}}+f{\mathit{e}}_z\wedge \mathit{u}+g{\mathbf{\nabla }}_{\mathrm{h}}\mathit{\eta }+{\mathit{F}}_{\mathrm{pg}}\\ & {\displaystyle }& {\displaystyle }={\mathbf{\nabla }}_{\mathrm{h}}\cdot \left({\mathit{\nu }}_{\mathrm{h}}{\mathbf{\nabla }}_{\mathrm{h}}\mathit{u}\right)+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\nu }{\displaystyle \frac{\partial \mathit{u}}{\partial z}}\right).\end{array}$$ Vertical velocity $w$ is diagnosed from the continuity equation: $${\displaystyle }{\displaystyle }{\mathbf{\nabla }}_{\mathrm{h}}\cdot \mathit{u}+{\displaystyle \frac{\partial w}{\partial z}}=\mathrm{0}.$$ Water temperature and salinity are modeled with an advection–diffusion equation of the form $${\displaystyle }{\displaystyle \frac{\partial T}{\partial t}}+{\mathbf{\nabla }}_{\mathrm{h}}\cdot \left(\mathit{u}T\right)+{\displaystyle \frac{\partial \left(wT\right)}{\partial z}}={\mathbf{\nabla }}_{\mathrm{h}}\cdot \left({\mathit{\mu }}_{\mathrm{h}}{\mathbf{\nabla }}_{\mathrm{h}}T\right)+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\mu }{\displaystyle \frac{\partial T}{\partial z}}\right),$$ where ${\mathit{\mu }}_{\mathrm{h}}$ and $\mathit{\mu }$ stand for the horizontal and vertical (eddy) diffusivity, respectively.

At the bottom boundary, we impose quadratic bottom stress: $$\begin{array}{}& {\displaystyle }& {\displaystyle }\left({\mathit{\nu }}_{\mathrm{h}}{\mathit{n}}_{\mathrm{h}}\cdot {\mathbf{\nabla }}_{\mathrm{h}}\mathit{u}+\mathit{\nu }{n}_z{\displaystyle \frac{\partial \mathit{u}}{\partial z}}\right)\left|{}_{x\in {\mathrm{\Gamma }}_{\mathrm{b}}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{b}}}{{\mathit{\rho }}_{\mathrm{0}}}},\right.& {\displaystyle }& {\displaystyle \frac{{\mathit{\tau }}_{\mathrm{b}}}{{\mathit{\rho }}_{\mathrm{0}}}}=C_{\mathrm{d}}\left|{\mathit{u}}_{\mathrm{bf}}\right|{\mathit{u}}_{\mathrm{bf}},\end{array}$$ where $C_{\mathrm{d}}$ is the drag coefficient, and ${\mathit{u}}_{\mathrm{bf}}$ is the velocity in the middle of the bottommost element. $\mathit{n}=(n_x$, $n_y$, $n_z)$ is the outward normal vector, and ${\mathit{n}}_{\mathrm{h}}=(n_x$, $n_y$, 0) its horizontal projection. The bottom boundary condition is treated implicitly; Eq. () is linearized by keeping the magnitude $\left|{\mathit{u}}_{\mathrm{bf}}\right|$ fixed at the “old” value while solving for $\mathit{u}$ (and ${\mathit{u}}_{\mathrm{bf}}$). Typically, $C_{\mathrm{d}}$ is computed from the logarithmic law of the wall e.g.,.

Mode splitting

Following , we split the horizontal velocity field into depth-averaged $\stackrel{\mathrm{‾}}{\mathit{u}}$ and deviation ${\mathit{u}}^{\prime }=\mathit{u}-\stackrel{\mathrm{‾}}{\mathit{u}}$ components. The depth-averaged momentum equation is then defined as $${\displaystyle }{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{u}}}{\partial t}}+f{\mathit{e}}_z\wedge \stackrel{\mathrm{‾}}{\mathit{u}}+g{\mathbf{\nabla }}_{\mathrm{h}}\mathit{\eta }=\mathit{G},$$ where $\mathit{G}$ is a forcing term used to couple the 2-D and 3-D modes. This equation is complemented with the depth-averaged continuity (free surface) equation: $${\displaystyle }{\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+{\mathbf{\nabla }}_{\mathrm{h}}\cdot \left(H\stackrel{\mathrm{‾}}{\mathit{u}}\right)=\mathrm{0}.$$ The 2-D system (Eqs. –) contains the fast-propagating, rotational surface gravity waves. The corresponding equation for ${\mathit{u}}^{\prime }$ is obtained by subtracting Eq. () from Eq. () : $$\begin{array}{cc}{\displaystyle \frac{\partial {\mathit{u}}^{\prime }}{\partial t}}& {\displaystyle }+{\mathbf{\nabla }}_{\mathrm{h}}\cdot \left(\mathit{u}\mathit{u}\right)+{\displaystyle \frac{\partial \left(w\mathit{u}\right)}{\partial z}}+f{\mathit{e}}_z\wedge {\mathit{u}}^{\prime }+{\mathit{F}}_{\mathrm{pg}}\\ & {\displaystyle }& {\displaystyle }={\mathbf{\nabla }}_{\mathrm{h}}\cdot \left({\mathit{\nu }}_{\mathrm{h}}{\mathbf{\nabla }}_{\mathrm{h}}\mathit{u}\right)+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\nu }{\displaystyle \frac{\partial \mathit{u}}{\partial z}}\right)-\mathit{G}.\end{array}$$ Note that the advection and viscosity terms are included in Eq. () without splitting, based on the assumption that these processes are slow enough to be captured with long time steps. The Coriolis term, on the other hand, only contains the slow modes. The vertical velocity $w$ only appears in the advection term, which is not split, and thus there is no need to split $w$.

Coupling 2-D and 3-D modes

The 2-D and 3-D modes are coupled using the additional term $\mathit{G}$ . First, the 3-D momentum equation (Eq. ) is solved with $\mathit{G}=\mathrm{0}$, resulting in a velocity field ${\mathit{u}}^{\prime }$ that has a non-zero depth average, generated by the advection and viscosity terms (that depend on $\stackrel{\mathrm{‾}}{\mathit{u}}$). We then compute the depth-averaged $\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }}$ and apply a correction: $$\begin{array}{}& {\displaystyle }& {\displaystyle }\mathit{G}=\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }}/\mathrm{\Delta }t,& {\displaystyle }& {\displaystyle }{\mathit{u}}^{\prime }\leftarrow {\mathit{u}}^{\prime }-\mathit{G}\mathrm{\Delta }t\end{array}$$ to enforce zero depth average. By definition, the field $\mathit{G}$ is a constant over the vertical, and it will be used as a forcing term in the 2-D momentum equation (Eq. ) in the subsequent solve. This procedure ensures that Eqs. () and () sum up to Eq. () and $\int {\mathit{u}}^{\prime }\mathrm{d}z\mathrm{0}$.

Equation of state

In this paper, a linear equation of state is used: $${\displaystyle }{\displaystyle }\mathit{\rho }(T,S)={\mathit{\rho }}_{\mathrm{0}}-{\mathit{\alpha }}_{\mathrm{T}}\left(T-T_{\mathrm{0}}\right)+{\mathit{\beta }}_{\mathrm{S}}\left(S-S_{\mathrm{0}}\right),$$ where ${\mathit{\alpha }}_{\mathrm{T}}$ and ${\mathit{\beta }}_{\mathrm{S}}$ are the thermal expansion and saline contraction coefficients, respectively, and $T_{\mathrm{0}}$ and $S_{\mathrm{0}}$ are reference temperature and salinity. In all the test cases presented herein, salinity does not contribute to water density (${\mathit{\beta }}_{\mathrm{S}}=\mathrm{0}$). Thetis also implements a full non-linear equation of state .

Viscosity and turbulence closure

Baroclinic flows require some form of viscosity to filter out grid-scale noise. In this paper, we only consider Laplacian horizontal viscosity, set to a constant ${\mathit{\nu }}_{\mathrm{h}}=U\mathrm{\Delta }x/$Re${}_{\mathrm{h}}$ corresponding to the velocity scale $U$, horizontal mesh resolution $\mathrm{\Delta }x$, and the desired grid Reynolds number Re${}_{\mathrm{h}}$. Here, the velocity scale $U$ is taken as a global constant specific to each test case. Unless otherwise specified, the horizontal diffusivity of tracers is zero.

In most test cases, vertical viscosity is set to a constant. In certain cases, we use the gradient Richardson number dependent parameterization by : $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathit{\nu }={\displaystyle \frac{{\mathit{\nu }}_{\mathrm{0}}}{(\mathrm{1}+\mathit{\alpha }\mathit{\text{Ri}}{)}^n}}+{\mathit{\nu }}_{\mathrm{b}},\\ & {\displaystyle }& {\displaystyle }\mathit{\mu }={\displaystyle \frac{\mathit{\nu }}{\mathrm{1}+\mathit{\alpha }\mathit{\text{Ri}}}}+{\mathit{\mu }}_{\mathrm{b}},\end{array}$$ where $\mathit{\text{Ri}}=N^{\mathrm{2}}/M^{\mathrm{2}}$ is the gradient Richardson number, $N$ is the buoyancy frequency, and $M$ is the vertical shear frequency. The background values are set to ${\mathit{\nu }}_{\mathrm{b}}={\mathit{\mu }}_{\mathrm{b}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$, while maximum viscosity is set to ${\mathit{\nu }}_{\mathrm{0}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{2}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$; the dimensionless parameters are $\mathit{\alpha }=\mathrm{10}$ and $n=\mathrm{2}$ . More sophisticated turbulence closures will be addressed in future work.

Finite element discretization

This section describes the spatial discretization of the governing equations. In Sect. , we define the finite element function spaces, followed by the weak forms of the underlying equations.

Prognostic and diagnostic variables and their function spaces.

Field	Symbol	Equation	Function space
Prognostic variables
Water elevation	$\mathit{\eta }$	()	$P_{\mathrm{1}}^{\mathrm{DG}}$
Depth av. velocity	$\stackrel{\mathrm{‾}}{\mathit{u}}$	()	$[P_{\mathrm{1}}^{\mathrm{DG}}{]}^{\mathrm{2}}$
Horizontal velocity	${\mathit{u}}^{\prime }$	()	$[P_{\mathrm{1}}^{\mathrm{DG}}\times P_{\mathrm{1}}^{\mathrm{DG}}{]}^{\mathrm{2}}$
Water temperature	$T$	()	$P_{\mathrm{1}}^{\mathrm{DG}}\times P_{\mathrm{1}}^{\mathrm{DG}}$
Water salinity	$S$	()	$P_{\mathrm{1}}^{\mathrm{DG}}\times P_{\mathrm{1}}^{\mathrm{DG}}$
Diagnostic variables
Vertical velocity	$w$	()	$P_{\mathrm{1}}^{\mathrm{DG}}\times P_{\mathrm{1}}^{\mathrm{DG}}$
Water density	${\mathit{\rho }}^{\prime }$	()	$P_{\mathrm{1}}^{\mathrm{DG}}\times P_{\mathrm{1}}^{\mathrm{DG}}$
Baroclinic head	$r$	()	$P_{\mathrm{1}}^{\mathrm{DG}}\times P_{\mathrm{2}}$
Int. pressure grad.	${\mathit{F}}_{\mathrm{pg}}$	()	$[P_{\mathrm{1}}^{\mathrm{DG}}\times P_{\mathrm{1}}^{\mathrm{DG}}{]}^{\mathrm{2}}$

Function spaces

The prognostic variables of the coupled 2-D–3-D system (Eqs. , , , ) are $\mathit{\eta }$, $\stackrel{\mathrm{‾}}{\mathit{u}}$, ${\mathit{u}}^{\prime }$, $T$, and $S$. Diagnostic variables include the vertical velocity $w$, water density ${\mathit{\rho }}^{\prime }$, baroclinic head $r$, and internal pressure gradient ${\mathit{F}}_{\mathrm{pg}}$. The choice of function spaces where these variables reside is crucial for numerical stability and accuracy.

Our discretization is based on the linear discontinuous Galerkin function space, $P_{\mathrm{1}}^{\mathrm{DG}}$. The 2-D system is discretized with a $P_{\mathrm{1}}^{\mathrm{DG}}-P_{\mathrm{1}}^{\mathrm{DG}}$ velocity–pressure finite element pair: water elevation and both components of the depth-averaged velocity are approximated in $P_{\mathrm{1}}^{\mathrm{DG}}$ space, i.e., $\mathit{\eta }\in {\mathcal{H}}_{\mathrm{2}\text{-}\mathrm{D}}=P_{\mathrm{1}}^{\mathrm{DG}}$, $\mathit{u}\in {\mathcal{U}}_{\mathrm{2}\text{-}\mathrm{D}}=[P_{\mathrm{1}}^{\mathrm{DG}}{]}^{\mathrm{2}}$. When embedded with appropriate Riemann fluxes at element interfaces, the $P_{\mathrm{1}}^{\mathrm{DG}}-P_{\mathrm{1}}^{\mathrm{DG}}$ element pair is well suited for rotational shallow water problems .

Achieving an accurate and monotone 3-D tracer advection scheme is one of our main design criteria. The tracers, therefore, are also considered within a discontinuous function space, $T$, $S\in \mathcal{H}=P_{\mathrm{1}}^{\mathrm{DG}}\times P_{\mathrm{1}}^{\mathrm{DG}}$ (here, the $\times$ operator stands for the Cartesian product of function spaces in the extruded mesh: horizontal $\times$ vertical function space). Tracer consistency (sometimes called local tracer conservation) is a necessary condition for monotonicity; it ensures that a constant tracer field does not exhibit spurious local extrema. In practice, it implies that the discrete tracer equation must reduce to the discrete continuity equation for a constant tracer. In this work, we satisfy this property by requiring the vertical velocity to belong to the tracer space $\mathcal{H}$ . In addition, compatibility between the 2-D and 3-D momentum equations requires that the 3-D horizontal velocity must be $P_{\mathrm{1}}^{\mathrm{DG}}$ in the horizontal direction. We therefore set ${\mathit{u}}^{\prime }\in \mathcal{U}=[P_{\mathrm{1}}^{\mathrm{DG}}\times P_{\mathrm{1}}^{\mathrm{DG}}{]}^{\mathrm{2}}$ as well. The used function spaces are summarized in Table .

Note that this choice of function spaces is not mimetic : the discrete system does not preserve all the properties of the continuous equations; for example, enstrophy is not conserved exactly. As the coastal ocean is generally very dissipative, maintaining mimetic properties is, however, not crucial. It is possible to define a mimetic discretization as well, for example, using Raviart–Thomas elements for the velocity, i.e., element pair RT${}_{\mathrm{1}}-P_{\mathrm{1}}^{\mathrm{DG}}$ . Our preliminary experiments, however, indicate that this choice significantly increases the computational cost of the system, without a corresponding improvement in accuracy. Formal assessment of the performance of mimetic discretizations in coastal ocean applications will be investigated in the future.

In the weak forms, we use the following notation for volume and interface integrals: $$\begin{array}{}& {\displaystyle }& {\displaystyle }⟨•{⟩}_{\mathrm{\Omega }}=\underset{\mathrm{\Omega }}{\int }•\mathrm{d}\mathbf{x},& {\displaystyle }& {\displaystyle }⟨⟨•⟩{⟩}_{\partial \mathrm{\Omega }}=\underset{\partial \mathrm{\Omega }}{\int }•\mathrm{d}s.\end{array}$$ In interface terms, we additionally use the average $\mathit{\left\{}\mathit{\right\{}\cdot \mathit{\left\}}\mathit{\right\}}$ and jump $\left[\right[\cdot \left]\right]$ operators for scalar $a$ and vector $\mathit{u}$ fields: $$\begin{array}{}& {\displaystyle }& {\displaystyle }\mathit{\left\{}\mathit{\right\{}a\mathit{\left\}}\mathit{\right\}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(a^{+}+a^{-}\right),& {\displaystyle }& {\displaystyle }\mathit{\left\{}\mathit{\right\{}\mathit{u}\mathit{\left\}}\mathit{\right\}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{u}}^{+}+{\mathit{u}}^{-}\right),& {\displaystyle }& {\displaystyle }\left[\right[a\mathit{n}\left]\right]=a^{+}{\mathit{n}}^{+}+a^{-}{\mathit{n}}^{-},& {\displaystyle }& {\displaystyle }\left[\right[\mathit{u}\cdot \mathit{n}\left]\right]={\mathit{u}}^{+}\cdot {\mathit{n}}^{+}+{\mathit{u}}^{-}\cdot {\mathit{n}}^{-},& {\displaystyle }& {\displaystyle }\left[\right[\mathit{u}\mathit{n}\left]\right]={\mathit{u}}^{+}{\mathit{n}}^{+}+{\mathit{u}}^{-}{\mathit{n}}^{-},\end{array}$$ where the superscripts “$+$” and “$-$” arbitrarily label the values on either side of the interface, and $\mathit{n}$ is the outward unit normal vector of each element on the interface.

2-D system

Let $\mathcal{T}$ stand for the triangulation of the 2-D domain ${\mathrm{\Gamma }}_{\mathrm{0}}$. The set of element interfaces is denoted by $\mathcal{I}=\mathit{\{}k\cap k^{\prime }\mathrm{|}k,k^{\prime }\in \mathcal{T}\mathit{\}}$, and $\mathit{n}=(n_x$, $n_y)$ the outward unit normal vector of an interface $e\in \mathcal{I}$. For brevity, boundary conditions are omitted from the weak forms.

Let ${\mathit{\varphi }}_{\mathrm{2}\text{-}\mathrm{D}}\in {\mathcal{H}}_{\mathrm{2}\text{-}\mathrm{D}}$ and ${\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}\in {\mathcal{U}}_{\mathrm{2}\text{-}\mathrm{D}}$ be test functions in the 2-D function spaces. The weak formulation of the 2-D system then reads: find $\mathit{\eta }\in {\mathcal{H}}_{\mathrm{2}\text{-}\mathrm{D}}$, $\stackrel{\mathrm{‾}}{\mathit{u}}\in {\mathcal{U}}_{\mathrm{2}\text{-}\mathrm{D}}$ such that $$\begin{array}{cc}{\displaystyle }⟨{\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}{\mathit{\varphi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}& {\displaystyle }+⟨⟨\left(H^{*}{\stackrel{\mathrm{‾}}{\mathit{u}}}^{*}\right)\cdot \left[\left[{\mathit{\varphi }}_{\mathrm{2}\text{-}\mathrm{D}}\mathit{n}\right]\right]⟩{⟩}_{\mathcal{I}}\\ & {\displaystyle }& {\displaystyle }-⟨\left(H\stackrel{\mathrm{‾}}{\mathit{u}}\right)\cdot {\mathbf{\nabla }}_{\mathrm{h}}{\mathit{\varphi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}=\mathrm{0},\end{array}$$ $$\begin{array}{cc}{\displaystyle }⟨{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{u}}}{\partial t}}\cdot {\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}& {\displaystyle }+⟨f{\mathit{e}}_z\wedge \stackrel{\mathrm{‾}}{\mathit{u}}\cdot {\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}+⟨⟨g{\mathit{\eta }}^{*}\left[\left[{\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}\cdot \mathit{n}\right]\right]⟩{⟩}_{\mathcal{I}}\\ {\displaystyle }& {\displaystyle }-⟨g\mathit{\eta }{\mathbf{\nabla }}_{\mathrm{h}}\cdot {\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}=⟨\mathit{G}\cdot {\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}},\\ & {\displaystyle }& {\displaystyle }\forall {\mathit{\varphi }}_{\mathrm{2}\text{-}\mathrm{D}}\in {\mathcal{H}}_{\mathrm{2}\text{-}\mathrm{D}},{\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}\in {\mathcal{U}}_{\mathrm{2}\text{-}\mathrm{D}}.\end{array}$$ Here, the divergence ${\mathbf{\nabla }}_{\mathrm{h}}\cdot \left(H\stackrel{\mathrm{‾}}{\mathit{u}}\right)$ and external gradient $g{\mathbf{\nabla }}_{\mathrm{h}}\mathit{\eta }$ terms have been integrated by parts. The resulting interface terms are defined on the element edges where the state variables $\mathit{\eta }$, $\stackrel{\mathrm{‾}}{\mathit{u}}$ are not uniquely defined. The values ${\mathit{\eta }}^{*}$, ${\stackrel{\mathrm{‾}}{\mathit{u}}}^{*}$ are obtained from an approximate Riemann solver; here, we use the linear Roe solution ${\mathit{\eta }}^{*}\mathit{\left\{}\mathit{\right\{}\mathit{\eta }\mathit{\left\}}\mathit{\right\}}+\sqrt{H/g}\left[\right[\stackrel{\mathrm{‾}}{\mathit{u}}\cdot \mathit{n}\left]\right]$ and ${\stackrel{\mathrm{‾}}{\mathit{u}}}^{*}=\mathit{\left\{}\mathit{\right\{}\stackrel{\mathrm{‾}}{\mathit{u}}\mathit{\left\}}\mathit{\right\}}+\sqrt{g/H}\left[\right[\mathit{\eta }\mathit{n}\left]\right]$ .

Momentum equation

Let $\mathcal{P}$ denote the set of prisms of the 3-D domain $\mathrm{\Omega }$, obtained from a vertical extrusion of ${\mathrm{\Gamma }}_{\mathrm{0}}$. The set of horizontal and vertical interfaces is denoted by ${\mathcal{I}}_{\mathrm{h}}$ and ${\mathcal{I}}_{\mathrm{v}}$, respectively. Let $\mathit{\psi }\in \mathcal{U}$ be a test function. The weak formulation of the 3-D momentum equation then reads: find $\mathit{u}\in \mathcal{U}$ such that $$\begin{array}{cc}{\displaystyle }⟨{\displaystyle \frac{\partial {\mathit{u}}^{\prime }}{\partial t}}\cdot \mathit{\psi }{⟩}_{\mathrm{\Omega }}& {\displaystyle }-⟨{\mathbf{\nabla }}_{\mathrm{h}}\mathit{\psi }:\left(\mathit{u}\mathit{u}\right){⟩}_{\mathrm{\Omega }}+⟨⟨{\mathit{u}}^{\mathrm{up}}\cdot \left[\left[\mathit{\psi }{\mathit{n}}_{\mathrm{h}}\right]\right]\cdot \mathit{\left\{}\mathit{\right\{}\mathit{u}\mathit{\left\}}\mathit{\right\}}⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}\cup {\mathcal{I}}_{\mathrm{v}}}\\ {\displaystyle }& {\displaystyle }-⟨\left(w\mathit{u}\right)\cdot {\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}{⟩}_{\mathrm{\Omega }}+⟨⟨{\mathit{u}}^{\mathrm{up}}\cdot \left[\left[\mathit{\psi }{n}_z\right]\right]\mathit{\left\{}\mathit{\right\{}w\mathit{\left\}}\mathit{\right\}}⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}}\\ {\displaystyle }& {\displaystyle }+⟨f{\mathit{e}}_z\wedge {\mathit{u}}^{\prime }\cdot \mathit{\psi }{⟩}_{\mathrm{\Omega }}+⟨{\mathit{F}}_{\mathrm{pg}}\cdot \mathit{\psi }{⟩}_{\mathrm{\Omega }}\\ {\displaystyle }& {\displaystyle }+⟨⟨{\mathit{\gamma }}_{\mathrm{lf}}\left[\right[\mathit{u}\left]\right]\cdot \left[\right[\mathit{\psi }\left]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}\cup {\mathcal{I}}_{\mathrm{v}}}=D_{\mathrm{h}}(\mathit{u},\mathit{\psi })\\ & {\displaystyle }& {\displaystyle }+D_{\mathrm{v}}(\mathit{u},\mathit{\psi }),\forall \mathit{\psi }\in \mathcal{U}.\end{array}$$ Here, the advection and viscosity terms have been integrated by parts see; the colon operator is the Frobenius inner product, $\mathbf{A}:\mathbf{B}=\sum _{i,j}{A}_{i,j}{B}_{i,j}$, and ${\mathit{u}}^{\mathrm{up}}$ stands for the upwind value at the interface. The internal pressure gradient term has been augmented with the Lax–Friedrichs flux with parameter ${\mathit{\gamma }}_{\mathrm{lf}}=\mathit{\left\{}\mathit{\right\{}\left|\mathit{u}\right|\mathit{\left\}}\mathit{\right\}}$. Adding such a flux is required to stabilize the internal pressure gradient: it reduces noise in the velocity field and decreases spurious numerical mixing in baroclinic applications. The $D_{\mathrm{h}}$, $D_{\mathrm{v}}$ terms denote the diffusion operators introduced later.

Tracer equation

The weak formulation of the tracer equations is derived analogously: find $T\in \mathcal{H}$ such that $$\begin{array}{cc}{\displaystyle }⟨{\displaystyle \frac{\partial T}{\partial t}}\mathit{\varphi }{⟩}_{\mathrm{\Omega }}& {\displaystyle }-⟨T\mathit{u}\cdot {\mathbf{\nabla }}_{\mathrm{h}}\mathit{\varphi }{⟩}_{\mathrm{\Omega }}+⟨⟨T^{\mathrm{up}}\left[\left[\mathit{\varphi }{\mathit{n}}_{\mathrm{h}}\right]\right]\cdot \mathit{\left\{}\mathit{\right\{}\mathit{u}\mathit{\left\}}\mathit{\right\}}⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}\cup {\mathcal{I}}_{\mathrm{v}}}\\ {\displaystyle }& {\displaystyle }-⟨\left(Tw\right){\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}{⟩}_{\mathrm{\Omega }}+⟨⟨T^{\mathrm{up}}\left[\left[\mathit{\varphi }{n}_z\right]\right]\mathit{\left\{}\mathit{\right\{}w\mathit{\left\}}\mathit{\right\}}⟩{⟩}_{{\mathcal{I}}_{\mathrm{v}}}\\ & {\displaystyle }& {\displaystyle }=D_{\mathrm{h}}(T,\mathit{\varphi })+D_{\mathrm{v}}(T,\mathit{\varphi }),\forall \mathit{\varphi }\in \mathcal{H}.\end{array}$$ Note that we do not employ the Lax–Friedrichs flux in the tracer equation.

Symmetric interior penalty stabilization

The presented discretization is unstable for elliptic operators, and the diffusion operators require additional stabilization. Here, we use the symmetric interior penalty Galerkin (SIPG) method . The SIPG formulation of the tracer diffusion operators read $$\begin{array}{cc}{\displaystyle }{D}_{\mathrm{h}}(T,\mathit{\varphi })& {\displaystyle }=-⟨{\mathit{\mu }}_{\mathrm{h}}\left({\mathbf{\nabla }}_{\mathrm{h}}\mathit{\varphi }\right)\cdot \left({\mathbf{\nabla }}_{\mathrm{h}}T\right){⟩}_{\mathrm{\Omega }}\\ {\displaystyle }& {\displaystyle }+⟨⟨\left\{\left\{{\mathit{\mu }}_{\mathrm{h}}{\mathbf{\nabla }}_{\mathrm{h}}T\right\}\right\}\cdot \left[\left[\mathit{\varphi }{\mathit{n}}_{\mathrm{h}}\right]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}\cup {\mathcal{I}}_{\mathrm{v}}}\\ {\displaystyle }& {\displaystyle }+⟨⟨\left\{\left\{{\mathit{\mu }}_{\mathrm{h}}{\mathbf{\nabla }}_{\mathrm{h}}\mathit{\varphi }\right\}\right\}\cdot \left[\left[T{\mathit{n}}_{\mathrm{h}}\right]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}\cup {\mathcal{I}}_{\mathrm{v}}}\\ & {\displaystyle }& {\displaystyle }-⟨⟨\mathit{\left\{}\mathit{\right\{}\mathit{\sigma }\mathit{\left\}}\mathit{\right\}}\left\{\left\{{\mathit{\mu }}_{\mathrm{h}}\right\}\right\}\left[\left[T{\mathit{n}}_{\mathrm{h}}\right]\right]\cdot \left[\left[\mathit{\varphi }{\mathit{n}}_{\mathrm{h}}\right]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}\cup {\mathcal{I}}_{\mathrm{v}}},\end{array}$$ $$\begin{array}{cc}{\displaystyle }{D}_{\mathrm{v}}(T,\mathit{\varphi })& {\displaystyle }=-⟨\mathit{\mu }{\displaystyle \frac{\partial T}{\partial z}}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}{⟩}_{\mathrm{\Omega }}+⟨⟨\left\{\left\{\mathit{\mu }{\displaystyle \frac{\partial T}{\partial z}}\right\}\right\}\left[\left[\mathit{\varphi }{n}_z\right]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}}\\ {\displaystyle }& {\displaystyle }+⟨⟨\left\{\left\{\mathit{\mu }{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}\right\}\right\}\left[\left[T{n}_z\right]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}}\\ & {\displaystyle }& {\displaystyle }-⟨⟨\mathit{\left\{}\mathit{\right\{}\mathit{\sigma }\mathit{\left\}}\mathit{\right\}}\mathit{\left\{}\mathit{\right\{}\mathit{\mu }\mathit{\left\}}\mathit{\right\}}\left[\left[T{n}_z\right]\right]\left[\left[\mathit{\varphi }{n}_z\right]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}}.\end{array}$$ For the viscosity terms, we get $$\begin{array}{cc}{\displaystyle }{D}_{\mathrm{h}}(\mathit{u},\mathit{\psi })& {\displaystyle }=-⟨{\mathit{\nu }}_{\mathrm{h}}\left({\mathbf{\nabla }}_{\mathrm{h}}\mathit{\psi }\right):{\left({\mathbf{\nabla }}_{\mathrm{h}}\mathit{u}\right)}^T{⟩}_{\mathrm{\Omega }}\\ {\displaystyle }& {\displaystyle }+⟨⟨\left[\left[\mathit{\psi }{\mathit{n}}_{\mathrm{h}}\right]\right]\cdot \left\{\left\{{\mathit{\nu }}_{\mathrm{h}}{\mathbf{\nabla }}_{\mathrm{h}}\mathit{u}\right\}\right\}⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}\cup {\mathcal{I}}_{\mathrm{v}}}\\ {\displaystyle }& {\displaystyle }+⟨⟨\left[\left[\mathit{u}{\mathit{n}}_{\mathrm{h}}\right]\right]\cdot \left\{\left\{{\mathit{\nu }}_{\mathrm{h}}{\mathbf{\nabla }}_{\mathrm{h}}\mathit{\psi }\right\}\right\}⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}\cup {\mathcal{I}}_{\mathrm{v}}}\\ & {\displaystyle }& {\displaystyle }-⟨⟨\mathit{\left\{}\mathit{\right\{}\mathit{\sigma }\mathit{\left\}}\mathit{\right\}}\left\{\left\{{\mathit{\nu }}_{\mathrm{h}}\right\}\right\}\left[\left[\mathit{u}{\mathit{n}}_{\mathrm{h}}\right]\right]\left[\left[\mathit{\psi }{\mathit{n}}_{\mathrm{h}}\right]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}\cup {\mathcal{I}}_{\mathrm{v}}},\end{array}$$ $$\begin{array}{cc}{\displaystyle }{D}_{\mathrm{v}}(\mathit{u},\mathit{\psi })& {\displaystyle }=-⟨\mathit{\nu }{\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}\cdot {\displaystyle \frac{\partial \mathit{u}}{\partial z}}{⟩}_{\mathrm{\Omega }}\\ {\displaystyle }& {\displaystyle }+⟨⟨\left[\left[\mathit{\psi }{n}_z\right]\right]\cdot \left\{\left\{\mathit{\nu }{\displaystyle \frac{\partial \mathit{u}}{\partial z}}\right\}\right\}⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}}\\ {\displaystyle }& {\displaystyle }+⟨⟨\left[\left[\mathit{u}{n}_z\right]\right]\cdot \left\{\left\{\mathit{\nu }{\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}\right\}\right\}⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}}\\ & {\displaystyle }& {\displaystyle }-⟨⟨\mathit{\left\{}\mathit{\right\{}\mathit{\sigma }\mathit{\left\}}\mathit{\right\}}\mathit{\left\{}\mathit{\right\{}\mathit{\nu }\mathit{\left\}}\mathit{\right\}}\left[\left[\mathit{u}{n}_z\right]\right]\cdot \left[\left[\mathit{\psi }{n}_z\right]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}}.\end{array}$$ The penalty factor $\mathit{\sigma }$ is defined as $\mathit{\sigma }=\mathit{\gamma }\frac{p(p+\mathrm{1})}{L}$ , where $p$ is the degree of the basis functions, $\mathit{\gamma }$ is a factor depending on mesh quality, and $L$ is the local element length scale in the normal direction of the interface. Let $h_{\mathrm{h}}$ and $h_{\mathrm{v}}$ denote the horizontal and vertical element sizes, and $\mathbf{\Delta }=\mathrm{diag}(h_{\mathrm{h}}$, $h_{\mathrm{h}}$, $h_{\mathrm{v}})$. We then define $L=\mathit{n}\cdot \mathbf{\Delta }\cdot \mathit{n}=h_{\mathrm{h}}(n_x^{\mathrm{2}}+n_y^{\mathrm{2}})+h_{\mathrm{v}}{n}_z^{\mathrm{2}}$ . In this paper, we use $\mathit{\gamma }=\mathrm{5}$.

Continuity equation

The vertical velocity $w$ is computed diagnostically from the continuity equation (Eq. ) by solving the weak form: find $w\in \mathcal{H}$ such that $$\begin{array}{cc}{\displaystyle }⟨⟨w{n}_z\mathit{\phi }⟩{⟩}_{{\mathrm{\Gamma }}_{\mathrm{s}}}& {\displaystyle }+⟨⟨\mathit{\left\{}\mathit{\right\{}w\mathit{\left\}}\mathit{\right\}}\left[\left[\mathit{\phi }{n}_z\right]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}}-⟨w{\displaystyle \frac{\partial \mathit{\phi }}{\partial z}}{⟩}_{\mathrm{\Omega }}\\ {\displaystyle }& {\displaystyle }=⟨\mathit{u}\cdot {\mathbf{\nabla }}_{\mathrm{h}}\mathit{\phi }{⟩}_{\mathrm{\Omega }}-⟨⟨\mathit{\left\{}\mathit{\right\{}\mathit{u}\mathit{\left\}}\mathit{\right\}}\cdot \left[\left[\mathit{\phi }{\mathit{n}}_{\mathrm{h}}\right]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}\cup {\mathcal{I}}_{\mathrm{v}}}\\ & {\displaystyle }& {\displaystyle }-⟨⟨\mathit{u}\cdot \mathit{\phi }{\mathit{n}}_{\mathrm{h}}⟩{⟩}_{{\mathrm{\Gamma }}_{\mathrm{s}}},\forall \mathit{\phi }\in \mathcal{H},\end{array}$$ where both the left- and right-hand sides have been integrated by parts. Note that the terms on the bottom surface ${\mathrm{\Gamma }}_{\mathrm{b}}$ vanish due to the impermeability constraint $\mathit{u}\cdot {\mathit{n}}_{\mathrm{h}}+w{n}_z=\mathrm{0}$.

Computing the internal pressure gradient

The water density is computed diagnostically using the equation of state. We use the same $P_{\mathrm{1}}^{\mathrm{DG}}\times P_{\mathrm{1}}^{\mathrm{DG}}$ function space for tracers and water density. In this work, we use a linear equation of state (Eq. ), and consequently density can be computed locally (at each node of the tracer field). In general, however, the equation of state is non-linear, and the density is projected on the $\mathit{\rho }$ field.

The baroclinic head is computed from Eq. () by integrating ${\mathit{\rho }}^{\prime }$ over the vertical. In practice, we solve equation $\frac{\partial r}{\partial z}={\mathit{\rho }}^{\prime }/{\mathit{\rho }}_{\mathrm{0}}$ weakly with the appropriate boundary conditions: $${\displaystyle }⟨⟨r{n}_z\mathit{\phi }⟩{⟩}_{{\mathrm{\Gamma }}_{\mathrm{s}}}+⟨⟨r_{\mathrm{up}}\left[\left[\mathit{\phi }{n}_z\right]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}}-⟨r{\displaystyle \frac{\partial \mathit{\phi }}{\partial z}}{⟩}_{\mathrm{\Omega }}=⟨{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathit{\rho }}^{\prime }\mathit{\phi }{⟩}_{\mathrm{\Omega }}.$$ Here, the left-hand side has been integrated by parts, and $r_{\mathrm{up}}$ denotes the value on the prism above the interface. Note that the free surface terms vanish because $r=\mathrm{0}$ on ${\mathrm{\Gamma }}_{\mathrm{s}}$ by definition. We use function space $P_{\mathrm{1}}^{\mathrm{DG}}\times P_{\mathrm{2}}$ for $r$ to alleviate internal pressure gradient errors .

Finally, taking a test function $\mathit{\psi }\in \mathcal{U}$, we compute the internal pressure gradient with the weak form $$\begin{array}{cc}{\displaystyle }⟨{\mathit{F}}_{\mathrm{pg}}\cdot \mathit{\psi }{⟩}_{\mathrm{\Omega }}& {\displaystyle }=-⟨gr{\mathbf{\nabla }}_{\mathrm{h}}\cdot \mathit{\psi }{⟩}_{\mathrm{\Omega }}+⟨⟨g\mathit{\left\{}\mathit{\right\{}r\mathit{\left\}}\mathit{\right\}}\left[\left[\mathit{\psi }\cdot {\mathit{n}}_{\mathrm{h}}\right]\right]⟩{⟩}_{{\mathcal{I}}_{\mathrm{h}}\cup {\mathcal{I}}_{\mathrm{v}}}\\ & {\displaystyle }& {\displaystyle }+⟨⟨gr\mathit{\psi }\cdot {\mathit{n}}_{\mathrm{h}}⟩{⟩}_{{\mathrm{\Gamma }}_{\mathrm{s}}\cup {\mathrm{\Gamma }}_{\mathrm{b}}},\forall \mathit{\psi }\in \mathcal{U},\end{array}$$ where the right-hand side has been integrated by parts. Usually, ${\mathit{F}}_{\mathrm{pg}}$ belongs to the same space as the horizontal velocity, i.e., $[P_{\mathrm{1}}^{\mathrm{DG}}\times P_{\mathrm{1}}^{\mathrm{DG}}{]}^{\mathrm{2}}$. However, to reduce bathymetry induced internal pressure gradient errors, it is possible to use a quadratic horizontal space, i.e., $r\in P_{\mathrm{2}}^{\mathrm{DG}}\times P_{\mathrm{2}}$ and ${\mathit{F}}_{\mathrm{pg}}\in [P_{\mathrm{2}}^{\mathrm{DG}}\times P_{\mathrm{1}}^{\mathrm{DG}}{]}^{\mathrm{2}}$. In this paper, we use a linear ${\mathit{F}}_{\mathrm{pg}}$ field unless otherwise specified.

Slope limiters

We use vertex-based $P_{\mathrm{1}}^{\mathrm{DG}}$ slope limiters for three-dimensional variables to ensure positivity. The limiter is applied to both tracer and horizontal velocity fields after each update of the advection operator as discussed in the next section.

Time integration

The coupled 2-D–3-D system is advanced in time with a two-stage ALE time integration scheme. In this section, we present the ALE formulation and summarize the final time integration scheme.

ALE mesh formulation

To accurately account for the free surface movement, one must move the mesh in the vertical direction. In this work, we adopt the ALE method . Here, we describe a mesh update procedure that stretches (or compresses) the mesh uniformly over the vertical direction. The ALE formulation, however, allows more complex mesh-moving methods as well, such as the (approximate) tracking of isopycnals .

In three dimensions, an ALE update consists of solving an advection–diffusion equation between two domains, ${\mathrm{\Omega }}^n$ and ${\mathrm{\Omega }}^{n+\mathrm{1}}$. Here, the domain is uniquely defined by the surface elevation field, such that for any time level $n$ the surface ${\mathrm{\Gamma }}_{\mathrm{s}}^n$ matches ${\mathit{\eta }}^n$. Due to the chosen discretization, the elevation field $\mathit{\eta }$ is discontinuous, yet we wish to maintain a conforming mesh, i.e., a continuous coordinate field $z$. This is achieved by projecting the elevation field ${\mathit{\eta }}^n$ to a continuous space and updating the geometry with the continuous field ${\mathit{\eta }}_{\mathrm{cg}}^n$. The projection induces a small discrepancy between the elevation field and the 3-D domain, but its effect remains negligible in practical applications because jumps in the elevation field are typically small.

Let ${\mathrm{\Omega }}^{\mathrm{ref}}$ be the reference domain corresponding to unperturbed elevation field ${\mathit{\eta }}_{\mathrm{cg}}=\mathrm{0}$, and $z_{\mathrm{ref}}\in [-h$, 0] its vertical coordinate. Applying a uniform mesh stretching, the time-dependent mesh coordinates can then be written as $${\displaystyle }{\displaystyle }{z}^n=z_{\mathrm{ref}}+{\mathit{\eta }}_{\mathrm{cg}}^n{\displaystyle \frac{z_{\mathrm{ref}}+h}{h}}\in \left[-h,{\mathit{\eta }}_{\mathrm{cg}}^n\right].$$ The mesh velocity is obtained as $w_{\mathrm{m}}=\frac{\partial z}{\partial t}$. In practice, the consecutive fields ${\mathit{\eta }}_{\mathrm{cg}}^{n+\mathrm{1}}$ and ${\mathit{\eta }}_{\mathrm{cg}}^n$ are known so we can evaluate $${\displaystyle }{\displaystyle }{w}_{\mathrm{m}}^{n+\mathrm{1}}={\displaystyle \frac{{\mathit{\eta }}_{\mathrm{cg}}^{n+\mathrm{1}}-{\mathit{\eta }}_{\mathrm{cg}}^n}{\mathrm{\Delta }t}}{\displaystyle \frac{z_{\mathrm{ref}}+h}{h}}.$$ Given the mesh velocity, a conservative ALE update can be written as $${\displaystyle }{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(⟨T\mathit{\varphi }{⟩}_{\mathrm{\Omega }}\right)=⟨F_{\mathrm{T}}\left(T,\mathit{u},w-w_{\mathrm{m}}\right)\mathit{\varphi }{⟩}_{\mathrm{\Omega }},$$ for a generic tracer equation $\frac{\partial T}{\partial t}=F_{\mathrm{T}}(T$, $\mathit{u}$, $w)$.

Coupled time integration scheme

The coupled 2-D–3-D system is advanced in time with a two-stage ALE time integration scheme. For convenience, we rewrite the 3-D momentum and tracer equations as $$\begin{array}{}& {\displaystyle }& {\displaystyle \frac{\partial T}{\partial t}}=F_{\mathrm{T}}(T,\mathit{u},w)+G_{\mathrm{T}}\left(T\right),& {\displaystyle }& {\displaystyle \frac{\partial \mathit{u}}{\partial t}}=F_{\mathit{u}}\left({\mathit{F}}_{\mathrm{pg}},\mathit{u},w\right)+G_{\mathit{u}}\left(\mathit{u}\right),\end{array}$$ where $F_{\mathrm{T}}$ and $F_{\mathit{u}}$ denote all the terms that are treated explicitly, while $G_{\mathrm{T}}$ and $G_{\mathit{u}}$ contain all the implicit terms. In this work, only vertical diffusion (Eq. ), vertical viscosity (Eq. ), and bottom friction terms are treated implicitly.

The explicit 3-D equations are advanced in time with a second-order SSP Runge–Kutta scheme, SSPRK(2,2) . For a generic problem ($\frac{\partial c}{\partial t}=F\left(c\right)$), the scheme reads $$\begin{array}{}& {\displaystyle }& {\displaystyle }{c}^{\left(\mathrm{1}\right)}=c^n+\mathrm{\Delta }tF\left(c^n\right),& {\displaystyle }& {\displaystyle }{c}^{n+\mathrm{1}}=c^n+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }tF\left(c^n\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }tF\left(c^{\left(\mathrm{1}\right)}\right).\end{array}$$

When applied to the explicit 3-D momentum and tracer equations, Eqs. () and (), both of these stages are ALE updates where the mesh is updated from geometry ${\mathrm{\Omega }}^n$ to ${\mathrm{\Omega }}^{\left(\mathrm{1}\right)}$ and then ${\mathrm{\Omega }}^{n+\mathrm{1}}$. The ALE formulation of the explicit 3-D tracer equation can then be written as $${\displaystyle }{\displaystyle }⟨T^{\left(\mathrm{1}\right)}\mathit{\varphi }{⟩}_{{\mathrm{\Omega }}^{\left(\mathrm{1}\right)}}=⟨T^n\mathit{\varphi }{⟩}_{{\mathrm{\Omega }}^n}+\mathrm{\Delta }t⟨F_{\mathrm{T}}\left(T^n,{\mathit{u}}^n,w^n-w_{\mathrm{m}}^{\left(\mathrm{1}\right)}\right)\mathit{\varphi }{⟩}_{{\mathrm{\Omega }}^n},$$ $$\begin{array}{cc}{\displaystyle }⟨{\stackrel{\mathrm{̃}}{T}}^{n+\mathrm{1}}\mathit{\varphi }{⟩}_{{\mathrm{\Omega }}^{n+\mathrm{1}}}& {\displaystyle }=⟨T^n\mathit{\varphi }{⟩}_{{\mathrm{\Omega }}^n}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }t⟨F_{\mathrm{T}}\left(T^n,{\mathit{u}}^n,w^n-w_{\mathrm{m}}^{\left(\mathrm{1}\right)}\right)\mathit{\varphi }{⟩}_{{\mathrm{\Omega }}^n}\\ & {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }t⟨F_{\mathrm{T}}\left(T^{\left(\mathrm{1}\right)},{\mathit{u}}^{\left(\mathrm{1}\right)},w^{\left(\mathrm{1}\right)}-w_{\mathrm{m}}^{n+\mathrm{1}}\right)\mathit{\varphi }{⟩}_{{\mathrm{\Omega }}^{\left(\mathrm{1}\right)}},\end{array}$$ where the vertical velocity is adjusted by the mesh velocity $w_{\mathrm{m}}$.

After the SSPRK update, the implicit terms are advanced with the backward Euler method. This step is computed in domain ${\mathrm{\Omega }}^{n+\mathrm{1}}$: $${\displaystyle }{\displaystyle }⟨T^{n+\mathrm{1}}\mathit{\varphi }{⟩}_{{\mathrm{\Omega }}^{n+\mathrm{1}}}=⟨{\stackrel{\mathrm{̃}}{T}}^{n+\mathrm{1}}\mathit{\varphi }{⟩}_{{\mathrm{\Omega }}^{n+\mathrm{1}}}+\mathrm{\Delta }t⟨G_{\mathrm{T}}\left(T^{n+\mathrm{1}}\right)\mathit{\varphi }{⟩}_{{\mathrm{\Omega }}^{n+\mathrm{1}}}.$$ The 3-D momentum equation is treated analogously.

The 2-D equations are advanced in time with an implicit scheme to circumvent the strict time step constraint imposed by surface gravity waves. To ensure consistency between the movement of the 3-D mesh and the 2-D mode, the 2-D time integration scheme must be compatible with the aforementioned SSPRK(2,2) method. Here, we use a combination of a forward Euler and trapezoidal steps: $$\begin{array}{}& {\displaystyle }& {\displaystyle }{c}^{\left(\mathrm{1}\right)}=c^n+\mathrm{\Delta }tF\left(c^n\right),& {\displaystyle }& {\displaystyle }{c}^{n+\mathrm{1}}=c^n+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }t\left(F\left(c^n\right)+F\left(c^{n+\mathrm{1}}\right)\right).\end{array}$$ Denoting the tendencies of the 2-D system (Eqs. –) by $F_{\mathit{\eta }}$ and $F_{\stackrel{\mathrm{‾}}{\mathit{u}}}$, respectively, we can write the 2-D solver as $$\begin{array}{}& {\displaystyle }& {\displaystyle }⟨{\mathit{\eta }}^{\left(\mathrm{1}\right)}{\mathit{\varphi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}=⟨{\mathit{\eta }}^n{\mathit{\varphi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}+\mathrm{\Delta }t⟨F_{\mathit{\eta }}\left({\mathit{\eta }}^n,{\stackrel{\mathrm{‾}}{\mathit{u}}}^n\right){\mathit{\varphi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}},& {\displaystyle }& {\displaystyle }⟨{\stackrel{\mathrm{‾}}{\mathit{u}}}^{\left(\mathrm{1}\right)}\cdot {\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}=⟨{\stackrel{\mathrm{‾}}{\mathit{u}}}^n\cdot {\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}+\mathrm{\Delta }t⟨F_{\stackrel{\mathrm{‾}}{\mathit{u}}}\left({\mathit{\eta }}^n,{\stackrel{\mathrm{‾}}{\mathit{u}}}^n\right)\cdot {\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}},\end{array}$$ $$\begin{array}{cc}{\displaystyle }⟨{\mathit{\eta }}^{n+\mathrm{1}}{\mathit{\varphi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}& {\displaystyle }=⟨{\mathit{\eta }}^n{\mathit{\varphi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}⟨\left(F_{\mathit{\eta }}\left({\mathit{\eta }}^n,{\stackrel{\mathrm{‾}}{\mathit{u}}}^n\right)\right.\\ & {\displaystyle }& {\displaystyle }\left..+F_{\mathit{\eta }}\left(H^n,{\stackrel{\mathrm{‾}}{\mathit{u}}}^{n+\mathrm{1}}\right)\right){\mathit{\varphi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{\mathrm{\Gamma }\mathrm{0}},\end{array}$$ $$\begin{array}{cc}{\displaystyle }⟨{\stackrel{\mathrm{‾}}{\mathit{u}}}^{n+\mathrm{1}}\cdot {\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}& {\displaystyle }=⟨{\stackrel{\mathrm{‾}}{\mathit{u}}}^n\cdot {\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}⟨\left(F_{\stackrel{\mathrm{‾}}{\mathit{u}}}\left({\mathit{\eta }}^n,{\stackrel{\mathrm{‾}}{\mathit{u}}}^n\right)\right.\\ & {\displaystyle }& {\displaystyle }\left..+F_{\stackrel{\mathrm{‾}}{\mathit{u}}}\left({\mathit{\eta }}^{n+\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathit{u}}}^{n+\mathrm{1}}\right)\right)\cdot {\mathit{\psi }}_{\mathrm{2}\text{-}\mathrm{D}}{⟩}_{{\mathrm{\Gamma }}_{\mathrm{0}}}.\end{array}$$ The second implicit stage is linearized by treating the total depth $H$ explicitly in Eq. ().

The 2-D system is solved first, resulting in an updated elevation field (${\mathit{\eta }}^{\left(\mathrm{1}\right)}$ and ${\mathit{\eta }}^{n+\mathrm{1}}$ for the two stages, respectively) and consequently mesh geometry (${\mathrm{\Omega }}^{\left(\mathrm{1}\right)}$ and ${\mathrm{\Omega }}^{n+\mathrm{1}}$). Once the mesh geometry is known, it is straightforward to compute the corresponding mesh velocity $w_{\mathrm{m}}$ and perform a 3-D ALE update.

The time integration method is second order for all the terms. The whole algorithm is summarized in Algorithm 1.

Choosing the time step

The maximal admissible time step is constrained by the stability of the coupled time integrator. The presented SSPRK(2,2) scheme has a CFL (Courant–Friedrichs–Lewy) factor 1. The 2-D scheme (Eq. ) and the implicit vertical solver (Eq. ), on the other hand, are unconditionally stable. This implies that the coupled system is stable under the same conditions as the explicit SSP scheme on its own.

 

[Figure omitted. See PDF]

The horizontal advection term imposes a constraint: $${\displaystyle }{\displaystyle }\mathrm{\Delta }{t}_{\mathrm{adv}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{h}}{L}_{\mathrm{h}}}{U}},$$ where $L_{\mathrm{h}}$ is the horizontal element size, $U$ is the maximal horizontal velocity scale, and ${\mathit{\sigma }}_{\mathrm{h}}$ is a length scaling factor. For the presented $P_{\mathrm{1}}^{\mathrm{DG}}$ discretization, we take $L_{\mathrm{h}}$ as the square root of the triangle area. For rectangular $P_{\mathrm{1}}^{\mathrm{DG}}$ elements and second-order RK schemes, the scaling factor is approximately ${\mathit{\sigma }}_{\mathrm{h}}=\mathrm{0.33}$ . In this work, we use ${\mathit{\sigma }}_{\mathrm{h}}=\mathrm{0.05}$ for all the diagnostic test cases. In strongly stratified flows, internal waves may impose a stricter constraint: $${\displaystyle }{\displaystyle }\mathrm{\Delta }{t}_{\mathrm{iw}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{h}}{L}_{\mathrm{h}}}{C_{\mathrm{iw}}+U}},$$ where $C_{\mathrm{iw}}$ is the speed of the internal waves.

Analogously, the time step constraint for vertical advection is $${\displaystyle }{\displaystyle }\mathrm{\Delta }{t}_{\mathrm{w}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{v}}{L}_z}{W}},$$ where $L_z$ is the element height, $W$ is the vertical velocity scale, and ${\mathit{\sigma }}_{\mathrm{v}}=\mathrm{0.125}$ is the scaling factor.

Given a horizontal viscosity scale $N_{\mathrm{h}}$, the explicit viscosity operator imposes a constraint: $${\displaystyle }{\displaystyle }\mathrm{\Delta }{t}_{\mathrm{visc}}={\mathit{\sigma }}_{\mathrm{visc}}{\displaystyle \frac{{\left({\mathit{\sigma }}_{\mathrm{h}}{L}_{\mathrm{h}}\right)}^{\mathrm{2}}}{N_{\mathrm{h}}}},$$ which may become stringent for small elements and large viscosity values. The scaling factor ${\mathit{\sigma }}_{\mathrm{visc}}$ depends on the used stabilization scheme; here, a value of ${\mathit{\sigma }}_{\mathrm{visc}}=\mathrm{2}$ is used. The constraint for horizontal diffusivity is analogous.

In the simulations presented herein, the minimal admissible time step is evaluated on the mesh based on constant a-priori velocity and viscosity scales. The time step is kept constant throughout the simulation.

Test cases

We demonstrate the performance of the proposed discretization with a suite of test cases of increasing complexity. We first evaluate the conservation and convergence of the solver in a barotropic standing wave test case. The convergence of baroclinic terms is then examined in a specific steady-state test case. The baroclinic solver and its numerical mixing are then evaluated with a non-rotating lock exchange test case and a rotating baroclinic eddy test, followed by the Dynamics of Overflow Mixing and Entrainment (DOME) overflow test.

Standing wave

We first evaluate the performance of the solver in a barotropic standing wave test case. The domain is a $L_x=\mathrm{60}$ km long rectangular channel, 625 m wide, and 100 m deep. All lateral boundaries are closed. Initially, the water is at rest. A 10 m tall sinusoidal elevation perturbation is prescribed along the channel (${\mathit{\eta }}_{\mathrm{a}}=-{\mathit{\eta }}_{\mathrm{0}}\cos (\mathrm{2}\mathit{\pi }x/L_x)$, ${\mathit{\eta }}_{\mathrm{0}}=\mathrm{10}$ m), leading to a non-linear wave as the simulation progresses.

The simulation is run for two full wave periods, approximately 3831.31 s. To investigate tracer conservation and consistency properties, two passive tracers are included: salinity is set to a constant 4.5 psu, while temperature varies between 5.0 and 15.0 ${}^{\circ }$C along the channel ($T=\mathrm{5}\sin (\mathrm{2}\mathit{\pi }x/L_x)+\mathrm{10}$ ${}^{\circ }$C).

The domain was discretized with a split-quad mesh using 40 elements along the channel (1500 m edge length) and four vertical layers. The time step is $\mathrm{\Delta }t=\mathrm{95.78}$ s, chosen to meet the horizontal advection condition.

During the simulation, the volume of the 3-D domain was conserved to accuracy $\mathcal{O}\left({\mathrm{10}}^{-\mathrm{15}}\right)$. The “2-D volume”, i.e., the integral of the elevation field, was conserved to accuracy $\mathcal{O}\left({\mathrm{10}}^{-\mathrm{16}}\right)$. Salinity remained at constant 4.5 psu with a small $\mathcal{O}\left({\mathrm{10}}^{-\mathrm{9}}\right)$ deviation. The total mass of salinity and temperature were both conserved to accuracy $\mathcal{O}\left({\mathrm{10}}^{-\mathrm{12}}\right)$. Over- and undershoots in the temperature field were negligible due to the slope limiters. Without the limiter, temperature overshoots were $\mathcal{O}\left({\mathrm{10}}^{-\mathrm{2}}\right)$. These results show that the model indeed fully conserves volume and tracers and does not exhibit overshoots. Moreover, the tracer consistency property is satisfied, verifying the integrity of the ALE formulation.

To investigate the order of convergence of the solver, we used a smaller initial elevation perturbation (${\mathit{\eta }}_{\mathrm{0}}=\mathrm{1}$ cm). In this case, the resulting standing wave is close to linear. At the end of the simulation, the solution was compared to the analytical solution of the linear wave equation (which coincides with the initial condition) by computing the $L_{\mathrm{2}}$ error, $\mathcal{E}\left(\mathit{\eta }\right)=\left(\underset{\mathrm{\Omega }}{\int }\right(\mathit{\eta }-{\mathit{\eta }}_{\mathrm{a}}{)}^{\mathrm{2}}\mathrm{d}\mathbf{x}{)}^{\mathrm{1}/\mathrm{2}}$.

Convergence of the $L_{\mathrm{2}}$ error in the standing wave test case. Tested element sizes were 3000, 1500, 1000, 750, 500, 375, and 300 m. The number indicates the slope of the least-squares best-fit line (dashed line).

[Figure omitted. See PDF]

We ran the simulation, varying the horizontal mesh resolution between 3 km and 300 m; the number of vertical levels varied between 2 and 20. In each case, the channel was made one element wide, and the time step was chosen to meet the CFL criterion for horizontal advection. At the end of the simulation, the $L_{\mathrm{2}}$ error was computed for water elevation and velocity (see Fig. ). The velocity field shows the expected second-order convergence, whereas elevation converges at a rate of 3.2. It is known that $P_{\mathrm{1}}^{\mathrm{DG}}$ shallow water equations models may exhibit superconvergence properties, especially for the elevation field . Here, our results verify that the solver behaves as expected and yields second-order accuracy under barotropic forcing.

Baroclinic MMS test

Verifying model accuracy under baroclinic forcing is more challenging as no analytical solutions are available. Here, we use the method of manufactured solutions MMS; to construct a steady-state test case that allows us to verify the correctness of the discrete baroclinic operators. The domain is a $L_x=\mathrm{15}$ by $L_y=\mathrm{10}$ km large and $h=\mathrm{40}$ m deep rectangular box. All lateral boundaries are closed. We prescribe initial velocity and temperature fields: $$\begin{array}{}& {\displaystyle }& {\displaystyle }{u}_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L_x}}x\right)\cos \left({\displaystyle \frac{\mathrm{3}z}{h}}\right),& {\displaystyle }& {\displaystyle }{v}_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\cos \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right)\sin \left({\displaystyle \frac{z}{\mathrm{2}h}}\right),& {\displaystyle }& {\displaystyle }{T}_{\mathrm{a}}=\mathrm{15}\sin \left({\displaystyle \frac{\mathit{\pi }x}{L_x}}\right)\sin \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right)\cos \left({\displaystyle \frac{z}{h}}\right)+\mathrm{15}.\end{array}$$ These functions were chosen to be analytic (infinitely differentiable) and fully three-dimensional to better quantify the spatial discretization error.

Salinity is set to a constant 35 psu. We use the linear equation of state (Eq. ) with ${\mathit{\rho }}_{\mathrm{0}}=\mathrm{1000}$ kg m${}^{-\mathrm{3}}$, ${\mathit{\alpha }}_{\mathrm{T}}=\mathrm{0.2}$ kg m${}^{-\mathrm{3}}$ ${}^{\circ }$C${}^{-\mathrm{1}}$, and $T_{\mathrm{0}}=\mathrm{5}$ ${}^{\circ }$C${}^{-\mathrm{1}}$. For the sake of simplicity, bathymetry is constant and elevation is set to zero initially. Coriolis frequency was set to a constant $f={\mathrm{10}}^{-\mathrm{4}}$ s${}^{-\mathrm{1}}$. Bottom friction, viscosity, and diffusivity are omitted.

Without any additional forcing, the initial conditions lead to a time-dependent solution. Following the MMS strategy, we add analytical source terms in the dynamic equations that cancel all the active terms in the equations, leading to a steady-state solution. The remaining error is purely the discretization error of the advection, pressure gradient, and Coriolis operators. The source terms are derived analytically and projected to the corresponding function space. The analytical formulae are given in Appendix .

The coarsest mesh contains four elements both in $x$ and $y$ directions and two vertical levels. We refine the mesh up to 10 times (40 elements and 20 vertical levels) and compute the $L_{\mathrm{2}}$ error of the prognostic fields against the exact solutions. In each case, the model is run for 50 iterations with a time step chosen to meet the CFL condition.

The variation of the $L_{\mathrm{2}}$ errors with resolution is shown in Fig. . All the prognostic variables exhibit the correct second-order convergence rate. The diagnostic vertical velocity field (which depends on the divergence of $\mathit{u}$) converges linearly as expected. Therefore, we conclude that advection, pressure gradient, and Coriolis terms are discretized correctly. We have also developed similar MMS tests for the diffusivity/viscosity operators and the bottom friction term, all of which show second-order convergence as well (not shown).

Convergence of the $L_{\mathrm{2}}$ error in the baroclinic MMS test case. The mesh was refined 1, 2, 4, 6, 8, and 10 times, resulting in resolutions of 2500, 1250, 625, 416.67, 312.5, and 250 m (shortest edge of the triangle). The time steps were 25.0, 12.5, 6.25, 4.167, 3.125, and 2.5 s, respectively. The number indicates the slope of the least-squares best-fit line (dashed line).

[Figure omitted. See PDF]

Lock exchange

The validity of the baroclinic solver and its level of spurious mixing is investigated with the standard lock exchange test case . Here, we follow the setup of and : The domain is a 64 km long and 1 km wide rectangular channel. Water depth is 20 m. Initially, the left-hand side of the domain is filled with dense water mass ($T=\mathrm{5.0}$ ${}^{\circ }$C) compared to the water on the right ($T=\mathrm{30.0}$ ${}^{\circ }$C). Salinity is kept at constant 35 psu. We use the same linear equation of state as in Sect. , resulting in a density difference of $\mathrm{\Delta }\mathit{\rho }=\mathrm{5.0}$ kg m${}^{-\mathrm{3}}$. The domain is discretized with a regular split-quad mesh. The triangle edge length is 500 m and 20 equidistant $\mathit{\sigma }$ levels are used in the vertical direction.

Stabilizing the internal pressure gradient requires some form of friction. To this end, we apply a constant Laplacian horizontal viscosity, using values $\mathit{\nu }=\mathrm{1.0}$, 10.0, 100.0, and 200.0 m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$. These values correspond to grid Reynolds numbers Re${}_{\mathrm{h}}=U\mathrm{\Delta }x/\mathit{\nu }=\mathrm{250.0}$, 25.0, 2.5, and 1.25, respectively, where the characteristic velocity scale of the exchange flow is $U=\mathrm{0.5}$ m s${}^{-\mathrm{1}}$. Vertical viscosity is set to a constant 10${}^{-\mathrm{4}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$. Bottom friction is omitted.

Figure  shows the initial density field and solution after 17 h of simulation for the three cases. Higher background viscosity leads to a less noisy velocity field and therefore sharper density front. The sharpness and shape of the fronts are similar to results presented in the literature e.g., Fig. 5 in. The low viscosity cases (Re${}_{\mathrm{h}}=\mathrm{25}\mathrm{250}$) exhibit an internal wave at the front which significantly increases the overall mixing.

Assuming that, in the absence of bottom friction, all available potential energy is transformed into kinetic energy, the maximum front propagation speed can be estimated as $c=\mathrm{1}/\mathrm{2}\sqrt{gH\mathrm{\Delta }\mathit{\rho }/{\mathit{\rho }}_{\mathrm{0}}}$ . Figure a shows the propagation of the front location at the bottom of the domain (the front at the surface behaves comparably). All three simulations are in good agreement with the theoretical propagation speed. The simulated front propagation speed is underestimated by roughly 5 %, indicating loss of energy due to mixing. These results are similar to results reported in the literature; e.g.,  show similar performance for ROMS, MITgcm, and MOM.

Figure b shows the maximum over- and undershoots in the temperature field during the simulation. Even in the low viscosity case (Re${}_{\mathrm{h}}=\mathrm{250}$), the overshoots are of order 10${}^{-\mathrm{5}}$ ${}^{\circ }$C, indicating that the tracer advection scheme is indeed close to monotone, due to the SSP time integration method and slope limiters. Note that, if the slope limiter is omitted, the overshoots can reach 30 ${}^{\circ }$C.

Water density in the lock exchange test case in the center of the domain ($y=\mathrm{0}$ km). (a) Initial condition. Solution after 17 h of simulation with Re${}_{\mathrm{h}}$ (b) 1.25, (c) 2.5, (d) 25.0, and (e) 250.0.

[Figure omitted. See PDF]

Diagnostics of the lock exchange test. (a) Location of the density front at the bottom of the domain, (b) over- and undershoots in the temperature field (with regard to to 30.0 and 5.0 ${}^{\circ }$C, respectively), and (c) normalized reference potential energy (RPE) versus simulation time.

[Figure omitted. See PDF]

To diagnose the role of spurious mixing, we use the reference potential energy (RPE; ). RPE is computed as the vertical center of mass of a sorted density field ${\mathit{\rho }}^{*}$: RPE $=$ $g\int {\mathit{\rho }}^{*}(z+h)\mathrm{d}\mathbf{x}$. The ${\mathit{\rho }}^{*}$ field is defined as the unique, stratified density field where the densest water parcels are distributed over the bottom, and density increases monotonically over the water column. As such, ${\mathit{\rho }}^{*}$ is the steady-state density distribution, and RPE represents the portion of potential energy that cannot be transformed into kinetic energy. Mixing the two water masses increases RPE (the center of mass), and thus the amount of unavailable potential energy increases. Figure c shows the evolution of normalized RPE, $\stackrel{\mathrm{‾}}{\mathrm{RPE}}\left(t\right)=\left(\mathrm{RPE}\right(t)-\mathrm{RPE}(\mathrm{0}\left)\right)/\mathrm{RPE}\left(\mathrm{0}\right)$, during the simulation. At $t=\mathrm{17}$ h, the values are 0.612, 1.13, 2.35, $\mathrm{3.11}\times {\mathrm{10}}^{-\mathrm{5}}$ for the four simulations. These results are in good agreement with those reported with MPAS-Ocean model : with the same mesh resolution, MPAS-Ocean shows slightly larger normalized RPE, for example, at $t=\mathrm{17}$ h $\stackrel{\mathrm{‾}}{\mathrm{RPE}}\approx \mathrm{3.5}\times {\mathrm{10}}^{-\mathrm{5}}$ in the case of Re${}_{\mathrm{h}}=\mathrm{25}$. The difference is likely due to the different spatial discretization ($P_{\mathrm{1}}^{\mathrm{DG}}$ instead of finite volumes) or differences in the numerical viscosity operator. Applying slope limiters to the velocity field is not necessary for numerical stability, but it reduces high-frequency noise in the velocity field and hence results in lower RPE values.

In order to investigate the role of the Lax–Friedrichs flux on numerical mixing, we ran the lock exchange test case with zero viscosity. After 17 h of simulation, the RPE value was approximately $\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{5}}$. When the Lax–Friedrichs flux was omitted, a similar RPE value was obtained with viscosity $\mathit{\nu }=\mathrm{3.125}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$. Therefore, in this particular test case, the Lax–Friedrichs flux introduces mixing that is roughly equivalent to 3 m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ viscosity, corresponding to Re${}_{\mathrm{h}}=\mathrm{80}$. When viscosity was non-zero, it was evident from the numerical simulations that the Lax–Friedrichs flux has a negligible impact on numerical mixing if Re $<$ 10 (not shown).

Baroclinic eddies

We investigate the model's ability to generate baroclinic eddies with the eddying channel test case of . This test case is an idealization of the Antarctic Circumpolar Current, the domain spanning 500 and 160 km in the meridional and zonal directions, respectively. The domain is 1000 m deep. At the zonal boundaries, periodic boundary conditions are applied; northern and southern boundaries are closed. The Coriolis parameter is taken as a constant $\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{4}}$ s${}^{-\mathrm{1}}$.

Initially, the domain is linearly stratified with warmer water at the surface. In addition, the northern half of the domain is warmer, with a narrow sinusoidal transition band separating the warm (northern) and cold (southern) water masses (Fig. ; see for the definition of the initial temperature field). Water temperature ranges between 10 and 20 ${}^{\circ }$C. A linear equation of state is used with ${\mathit{\rho }}_{\mathrm{0}}=\mathrm{1000}$ kg m${}^{-\mathrm{3}}$, ${\mathit{\alpha }}_{\mathrm{T}}=\mathrm{0.2}$ kg m${}^{-\mathrm{3}}$ ${}^{\circ }$C${}^{-\mathrm{1}}$ and $T_{\mathrm{0}}=\mathrm{5}$ ${}^{\circ }$C. Salinity is kept at constant 35 psu and it does not affect density (${\mathit{\beta }}_{\mathrm{S}}=\mathrm{0}$). Bottom friction is parameterized by a constant drag coefficient of $C_{\mathrm{D}}=\mathrm{0.01}$.

The baroclinic Rossby radius of deformation is 20 km . Horizontal mesh resolution is constant in space. We use a regular split-quad mesh with two different mesh resolutions: eddy-permitting 10 km and a finer 4 km resolution. In the vertical direction, 26 and 40 equidistant sigma levels are used in the two cases, respectively. Simulations are carried out with different values of horizontal viscosity, with the grid Reynolds number ranging from 2 to 100. The different setups are summarized in Table . Vertical viscosity is set to a constant 10${}^{-\mathrm{4}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$.

As the simulation progresses, baroclinic eddies develop at the center of the domain, quickly propagating elsewhere. This is a spin-down experiment, i.e., the domain is a closed system with no forcing at the boundaries. Therefore, all the energy in the system originates from the initial potential energy, which is being dissipated during the simulation; again, the RPE is used as a metric for the energy transfer or the loss of energy due to mixing.

Experimental setup for baroclinic eddy test case. Listed are the horizontal mesh resolution (min. triangle edge length), number of vertical levels, time step, horizontal viscosity, and the approximate grid Reynolds number.

$\mathrm{\Delta }x$	$nz$	$\mathrm{\Delta }t$	${\mathit{\nu }}_{\mathrm{h}}$	Re${}_{\mathrm{h}}$
(km)		(s)	(m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$)
10	26	348.39	10.0	100
10	26	348.39	20.0	50
10	26	348.39	50.0	20
10	26	348.39	125.0	8
10	26	348.39	200.0	5
10	26	348.39	500.0	2
4	40	140.26	4.0	100
4	40	140.26	8.0	50
4	40	140.26	20.0	20
4	40	140.26	50.0	8
4	40	140.26	200.0	2

Figure  shows the surface temperature fields at various time intervals up to 200 days after the initialization for different values of horizontal viscosity. As expected, the model captures more mesoscale features as viscosity is decreased. Qualitatively, the results are in agreement with ROMS and MITgcm results , as well as MPAS-Ocean , all of which use a comparable Laplacian scheme for horizontal viscosity.

The evolution of the normalized RPE during the simulation is shown in Fig. a for the 4 km mesh resolution. The amount of mixing clearly depends on the grid Reynolds number, with RPE being roughly twice as high for Re${}_{\mathrm{h}}=\mathrm{20}$ compared to Re${}_{\mathrm{h}}=\mathrm{2}$. The average rate of change of RPE, averaged over days 3 to 319, is shown in Fig. b for all the simulations. As expected, the rate of change increases with larger grid Reynolds number and with a coarser mesh. These RPE metrics are in good agreement with results in the literature. At Re${}_{\mathrm{h}}=\mathrm{20}$ Thetis $\mathrm{dRPE}/\mathrm{d}t$, values are $\mathrm{4.3}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{4}}$ W m${}^{-\mathrm{2}}$, for the 10 and 4 km resolutions. The corresponding values for MITgcm, Modular Ocean Model (MOM), and Parallel Ocean Program (POP) (averaged over days 3 to 319) are larger: at least $\mathrm{8}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ W m${}^{-\mathrm{2}}$, respectively Fig. 12. reported similar values for MITgcm and MOM. On a hexagonal mesh, MPAS-Ocean yields smaller $\mathrm{dRPE}/\mathrm{d}t$ values: approximately $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{7}\times {\mathrm{10}}^{-\mathrm{5}}$ W m${}^{-\mathrm{2}}$ for the two resolutions, respectively values averaged over days 1–320; see Fig. 12 in. With a quad mesh, however, MPAS-Ocean values are approximately $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ W m${}^{-\mathrm{2}}$ for both resolutions and therefore close to Thetis performance.

Sea surface temperature fields for the eddying channel test case at 4 km horizontal mesh resolution. Horizontal viscosity is 200 (a), 50 (b), and 20 m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ (c). These values correspond to mesh Reynolds numbers 2, 8, and 20, respectively.

[Figure omitted. See PDF]

The test cases were run on a Linux cluster with 16-core Intel Xeon E5620 processors and Mellanox Infiniband interconnect. The 320-day simulation took roughly 42 h to run on 96 cores with the 4 km resolution mesh and 140.26 s time step. It should be noted, however, that the time step employed here is smaller than the maximal allowed time step. We also carried out a strong scaling test with the 4 km mesh. In the scaling test, the simulation was run for 40 time steps, recording the total elapsed wall-clock time and time spent in different parts of the solver. Figure a shows the overall speed-up up to 96 processors. The scaling efficiency drops to roughly 50 % at 96 cores, when the local degree of freedom count for the tracer field is 25 000 (see black line in Fig. b). This scaling efficiency is close to typical Firedrake performance .

Diagnostics of the eddying channel test case. (a) Evolution of normalized RPE over time in the eddying channel test case for 4 km mesh resolution. (b) Rate of change of RPE for different grid resolutions and grid Reynolds numbers. The rate of change was evaluated by computing the average RPE change from days 3 to 319.

[Figure omitted. See PDF]

Strong parallel scaling for the baroclinic eddies test case on a 4 km mesh (${\mathit{\nu }}_{\mathrm{h}}=\mathrm{20}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$): (a) speed-up in wall-clock time versus number of processes; (b) parallel efficiency versus the local number of degrees of freedom (DOFs) in the 3-D tracer field (top axis) and the 2-D ($\stackrel{\mathrm{‾}}{\mathit{u}}$, $\mathit{\eta }$) mixed system (bottom axis). The black line is the wall-clock time; colored lines stand for the time spent in different implicit or explicit solvers. The vertical dashed lines indicate 20 000 DOFs per process for the 2-D and 3-D problems, respectively. The mesh consisted of 10 000 triangles, 40 vertical levels, and 400 000 prisms.

[Figure omitted. See PDF]

The scaling efficiency of the separate solvers is plotted with colored lines in Fig. b. The implicit vertical diffusion/viscosity solvers perform best due to the fact that the problem is purely local without any horizontal dependencies. The explicit momentum solver scales almost as well, whereas the explicit tracer solver scales worse. The implicit 2-D solver (assembly and linear solve) scales the poorest because the problem is relatively small; at 96 cores, there are only around 940 degrees of freedom in the ($\stackrel{\mathrm{‾}}{\mathit{u}}$, $\mathit{\eta }$) system per core. We have also experimented with explicit 2-D solvers, but they do not scale significantly better compared to the two-stage implicit scheme used herein.

To further assess the CPU cost, we compared Thetis timing against the SLIM 3-D model which uses a similar DG formulation but is implemented in C/C$++$. The wall-clock time, and parallel efficiency used by both Thetis and SLIM 3-D are presented in Appendix . The setup, mesh, and time step were identical for the two models. On a single core, Thetis runs 3.3 times faster. On 24 cores, the ratio is 4.0, and on 144 cores Thetis is still 2.2 times faster than SLIM 3-D. This highlights the fact that Firedrake can deliver good parallel performance compared to models written in lower level languages.

It should be noted, however, that Thetis performance is currently not fully optimized. We expect that the performance can be significantly improved both in terms of serial and strong scaling performance. These will be addressed in future work.

Horizontal mesh and bathymetry for the DOME test case. The domain is extended 120 km further to the west to avoid boundary effects (shaded region). Horizontal element size ranges from 6 to 22 km. There are $\mathrm{18.8}\times {\mathrm{10}}^{\mathrm{3}}$ triangles in the horizontal mesh and 24 uniformly distributed vertical levels resulting in $\mathrm{450}\times {\mathrm{10}}^{\mathrm{3}}$ prisms and $\mathrm{2.7}\times {\mathrm{10}}^{\mathrm{6}}$ tracer degrees of freedom.

[Figure omitted. See PDF]

DOME

Next, we investigate the model's ability to simulate density-driven overflows with the DOME benchmark . The domain is a 1100 by 600 km large basin, whose depth varies linearly from 600 m at the northern boundary to 3600 m in the middle of the domain (see Fig. ). To avoid boundary condition issues, we have extended the domain to the west by 120 km. At the northern boundary, there is a 100 km wide and 200 km long inlet. Initially, the basin is stably stratified with a linear temperature variation from 10 ${}^{\circ }$C in the deepest part of the basin to 20 ${}^{\circ }$C at the surface. We use the linear equation of state with ${\mathit{\rho }}_{\mathrm{0}}=\mathrm{1000}$ kg m${}^{-\mathrm{3}}$, ${\mathit{\alpha }}_{\mathrm{T}}=\mathrm{0.2}$ kg m${}^{-\mathrm{3}}$ ${}^{\circ }$C${}^{-\mathrm{1}}$, and $T_{\mathrm{0}}=\mathrm{10}$ ${}^{\circ }$C, resulting in a $\mathrm{\Delta }\mathit{\rho }=\mathrm{2.0}$ kg m${}^{-\mathrm{3}}$ density difference.

At the inlet, a dense inflow (temperature 10 ${}^{\circ }$C) is prescribed in the bottom layer, with the surface layer being at 20 ${}^{\circ }$C. The inflow is in geostrophic balance, the thickness of the bottom layer being roughly 300 m on the eastern end of the boundary diminishing exponentially westward . The total inflow in the bottom layer is 5 Sv ($\mathrm{5}\times {\mathrm{10}}^{\mathrm{6}}$ m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$), the surface layer being static. During the simulation, the fate of the inflowing waters is tracked with a passive tracer that is initially zero in the basin and unity at the inlet. Initially, the tracer field is set to the inflow conditions in the northern part of the basin ($y>\mathrm{650}$ km). Velocity is set to zero everywhere. The eastern and southern boundaries of the basin are closed. The western boundary is open with radiation boundary conditions and a 100 km wide band where the temperature is relaxed to the initial condition.

The domain is discretized with an unstructured grid (Fig. ). Horizontal mesh resolution is 6 km near the northern boundary, increasing southward. Overall, 24 vertical sigma levels are used. Over the slope, the mesh resolution was designed to result in a hydrostatic consistency metric ($r<\mathrm{1.5}$) . Horizontal viscosity is set to a constant 50 m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$, which corresponds to Re${}_{\mathrm{h}}\approx \mathrm{200}$ at the inlet. Horizontal diffusivity is constant at 10 m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$. Vertical viscosity and diffusivity are parameterized by the Pacanowski–Philander scheme as described in Sect. . Bottom friction is parameterized with a quadratic drag coefficient $C_{\mathrm{d}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ . A quadratic function space is used for the baroclinic head and internal pressure gradient as discussed in Sect. .

As the inflowing current reaches the basin, it turns to the west and forms a coastal plume that is approximately 150 km wide (Fig. ). The plume detaches from the lateral boundary as it flows westward and along the bottom slope. As the dense water mass meets the stratified ocean, the plume becomes unstable and starts to generate eddies and internal waves. The most vigorous eddies are found in the first 300 km after the inlet ($x=\mathrm{500}$–800 km), after which the plume is more mixed and quiescent. Overall, the plume is shallow; most of the passive tracer is concentrated within 200 m of the bottom. Qualitatively, the extent and propagation of the plume, and its eddy structure are in good agreement with the literature e.g.,. The results show that Thetis is able to represent eddying flows over sloping bathymetry, generating and maintaining strong gradients between water masses. The sharpest fronts in the simulation encompass only one or two elements.

Figure  shows the distribution of the inflowing tracer concentration as a function of water density and the $x$ axis. The inflowing waters are initially very dense but get mixed to lower density as the plume advances along the coast. The histogram shows that the plume volume is low in the first 150 km after the inlet ($x=\mathrm{650}$–800 km) where the plume accelerates. After $x=\mathrm{650}$ km, the plume slows down and starts to accumulate in volume. The density of the main plume occupies ranges from 0.8 to 1.5 kg m${}^{-\mathrm{3}}$, the peak being around 1.28 kg m${}^{-\mathrm{3}}$. The rate of entrainment can be used as a metric for mixing. Results herein are similar to those presented in literature: present a mean density anomaly of 1.5 kg m${}^{-\mathrm{3}}$ for their terrain-following FESOM model configuration.

The 47-day simulation took roughly 42 h to run on 90 cores with a 39.65 s time step on the same Linux cluster.

Bottom tracer concentration in the DOME test case after 10 (a), 20 (b), and 40 days (c).

[Figure omitted. See PDF]

Histogram of tracer in the DOME test case versus the $x$ coordinate and density class. At the mouth of the inlet ($x=\mathrm{800}$ km), the inflowing waters are dense; they get entrained higher up in the density spectrum as they are being transported downstream. The data are averaged over 1 week after day 40.

[Figure omitted. See PDF]

Conclusions

This paper describes a DG implementation of an eddy-permitting, unstructured grid coastal ocean model. The solver is second-order accurate in space and time. We have demonstrated that the formulation is fully conservative and preserves monotonicity. The test cases indicate that the model is capable of reproducing the expected physical behavior, including baroclinic eddies. Moreover, numerical mixing is well-controlled and comparable to other established structured grid models, such as MITgcm and ROMS, and the large-scale finite volume model MPAS-Ocean. Finding an accurate formulation is important, as commonly used unstructured grid models tend to be overly diffusive, preventing accurate modeling of certain coastal domains e.g.,. The formulation presented herein thus contributes to the development of more accurate next-generation coastal ocean models.

Future work will include solving the equations on a sphere, DG implementation of a biharmonic viscosity operator, two-equation turbulence closure models, wetting–drying treatment, and development of an adjoint solver, as well as improving the computational efficiency and parallel scaling of the solver.

All code used to perform the experiments in this papers is publicly available. Firedrake, and its components, may be obtained from https://www.firedrakeproject.org/ (last access: 25 October 2018); Thetis from http://thetisproject.org/ (last access: 25 October 2018).

For reproducibility, we also cite archives of the exact software versions used to produce the results in this paper. All major Firedrake components have been archived on Zenodo . This record collates DOIs for the components and can be installed following the instructions at https://www.firedrakeproject.org/download.html (last access: 25 October 2018). Thetis itself has been archived at .

No external data were used in this paper.

Source terms for the baroclinic MMS test

Using the analytical velocity and temperature fields, we can derive the steady-state solution for the remaining fields: $$\begin{array}{}& {\displaystyle }& {\displaystyle }{\mathit{\eta }}_{\mathrm{a}}=\mathrm{0},& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{u}}_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}\sin \left(\mathrm{3}\right)\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L_x}}x\right),& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{v}}_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\sin \left({\displaystyle \frac{z}{\mathrm{2}h}}\right)\cos \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right),& {\displaystyle }& {\displaystyle }{u}_{\mathrm{a}}^{\prime }=u_{\mathrm{a}}-{\stackrel{\mathrm{‾}}{u}}_{\mathrm{a}},& {\displaystyle }& {\displaystyle }{v}_{\mathrm{a}}^{\prime }=v_{\mathrm{a}}-{\stackrel{\mathrm{‾}}{v}}_{\mathrm{a}},\end{array}$$ $$\begin{array}{cc}{\displaystyle }{w}_{\mathrm{a}}& {\displaystyle }={\displaystyle \frac{\mathit{\pi }h}{\mathrm{3}{L}_x{L}_y}}\left(\mathrm{2}{L}_x\left(-\cos \left({\displaystyle \frac{z}{\mathrm{2}h}}\right)+\cos \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)\right)\sin \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right)\right.\\ & {\displaystyle }& {\displaystyle }\left..-L_y\left(\sin \left({\displaystyle \frac{\mathrm{3}z}{h}}\right)+\sin \left(\mathrm{3}\right)\right)\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L_x}}x\right)\right),\end{array}$$ $${\displaystyle }{\displaystyle }{r}_{\mathrm{a}}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{T}}}{{\mathit{\rho }}_{\mathrm{0}}}}\left(T_{\mathrm{0}}z-\mathrm{15}h\sin \left({\displaystyle \frac{z}{h}}\right)\sin \left({\displaystyle \frac{\mathit{\pi }x}{L_x}}\right)\sin \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right)-\mathrm{15}z\right).$$

Now, we can evaluate the different terms that appear in the momentum and tracer equations: $$\begin{array}{}& {\displaystyle }& {\displaystyle }f{\left({\mathit{e}}_z\wedge \stackrel{\mathrm{‾}}{\mathit{u}}\right)}_x={\displaystyle \frac{\mathrm{2}{f}_{\mathrm{0}}}{\mathrm{3}}}\left(-\cos \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)+\mathrm{1}\right)\cos \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right),& {\displaystyle }& {\displaystyle }f{\left({\mathit{e}}_z\wedge \stackrel{\mathrm{‾}}{\mathit{u}}\right)}_y={\displaystyle \frac{f_{\mathrm{0}}}{\mathrm{6}}}\sin \left(\mathrm{3}\right)\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L_x}}x\right),\end{array}$$ $$\begin{array}{cc}{\displaystyle }{\mathbf{\nabla }}_{\mathrm{h}}\cdot \left(H\stackrel{\mathrm{‾}}{\mathit{u}}\right)& {\displaystyle }={\displaystyle \frac{\mathit{\pi }h}{\mathrm{3}{L}_x{L}_y}}\left(\mathrm{2}{L}_x\left(-\cos \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)\right.\right.\\ & {\displaystyle }& {\displaystyle }\left..\left..+\mathrm{1}\right)\sin \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right)+L_y\sin \left(\mathrm{3}\right)\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L_x}}x\right)\right),\end{array}$$ $$\begin{array}{cc}& {\displaystyle }& {\displaystyle }{\left({\mathit{F}}_{\mathrm{pg}}\right)}_x={\displaystyle \frac{\mathrm{15}\mathit{\pi }{\mathit{\alpha }}_{\mathrm{T}}h}{L_x{\mathit{\rho }}_{\mathrm{0}}}}g\sin \left({\displaystyle \frac{z}{h}}\right)\sin \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right)\cos \left({\displaystyle \frac{\mathit{\pi }x}{L_x}}\right),& {\displaystyle }& {\displaystyle }{\left({\mathit{F}}_{\mathrm{pg}}\right)}_y={\displaystyle \frac{\mathrm{15}\mathit{\pi }{\mathit{\alpha }}_{\mathrm{T}}h}{L_y{\mathit{\rho }}_{\mathrm{0}}}}g\sin \left({\displaystyle \frac{z}{h}}\right)\sin \left({\displaystyle \frac{\mathit{\pi }x}{L_x}}\right)\cos \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right),& {\displaystyle }& {\displaystyle }{\left({\mathbf{\nabla }}_{\mathrm{h}}\cdot \left(\mathit{u}\mathit{u}\right)\right)}_x={\displaystyle \frac{\mathit{\pi }}{\mathrm{2}{L}_x}}\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L_x}}x\right){\cos }^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{3}z}{h}}\right)\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L_x}}x\right),{\displaystyle }& {\displaystyle }{\left({\mathbf{\nabla }}_{\mathrm{h}}\cdot \left(\mathit{u}\mathit{u}\right)\right)}_y=-{\displaystyle \frac{\mathit{\pi }}{\mathrm{9}{L}_y}}{\sin }^{\mathrm{2}}\left({\displaystyle \frac{z}{\mathrm{2}h}}\right)\\ & {\displaystyle }& {\displaystyle }\sin \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right)\cos \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right),\end{array}$$ $$\begin{array}{cc}{\displaystyle \frac{\partial \left(wu\right)}{\partial z}}& {\displaystyle }={\displaystyle \frac{\mathit{\pi }\sin \left(\frac{\mathrm{3}z}{h}\right)}{\mathrm{2}{L}_x{L}_y}}\left(\mathrm{2}{L}_x\left(\cos \left({\displaystyle \frac{z}{\mathrm{2}h}}\right)-\cos \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)\right)\right.\\ {\displaystyle }& {\displaystyle }.\sin \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right)+L_y\left(\sin \left({\displaystyle \frac{\mathrm{3}z}{h}}\right)+\sin \left(\mathrm{3}\right)\right)\\ & {\displaystyle }& {\displaystyle }\left..\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L_x}}x\right)\right)\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L_x}}x\right),\end{array}$$ $$\begin{array}{cc}{\displaystyle \frac{\partial \left(wv\right)}{\partial z}}& {\displaystyle }=-{\displaystyle \frac{\mathit{\pi }\cos \left(\frac{z}{\mathrm{2}h}\right)}{\mathrm{18}{L}_x{L}_y}}\left(\mathrm{2}{L}_x\left(\cos \left({\displaystyle \frac{z}{\mathrm{2}h}}\right)-\cos \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)\right)\right.\\ {\displaystyle }& {\displaystyle }.\sin \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right)+L_y\left(\sin \left({\displaystyle \frac{\mathrm{3}z}{h}}\right)+\sin \left(\mathrm{3}\right)\right)\\ & {\displaystyle }& {\displaystyle }\left..\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L_x}}x\right)\right)\cos \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right),\end{array}$$ $$\begin{array}{cc}{\displaystyle }f{\left({\mathit{e}}_z\wedge {\mathit{u}}^{\prime }\right)}_x=& {\displaystyle \frac{f_{\mathrm{0}}}{\mathrm{3}}}\left(-\sin \left({\displaystyle \frac{z}{\mathrm{2}h}}\right)-\mathrm{2}+\mathrm{2}\cos \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)\right)\\ & {\displaystyle }& {\displaystyle }\cos \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right),\end{array}$$ $${\displaystyle }{\displaystyle }f{\left({\mathit{e}}_z\wedge {\mathit{u}}^{\prime }\right)}_y={\displaystyle \frac{f_{\mathrm{0}}}{\mathrm{6}}}\left(\mathrm{3}\cos \left({\displaystyle \frac{\mathrm{3}z}{h}}\right)-\sin \left(\mathrm{3}\right)\right)\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L_x}}x\right),$$ $$\begin{array}{cc}{\displaystyle }{\mathbf{\nabla }}_{\mathrm{h}}\cdot \left(\mathit{u}T\right)& {\displaystyle }={\displaystyle \frac{\mathrm{5}\mathit{\pi }}{L_x{L}_y}}\left(L_x\sin \left({\displaystyle \frac{z}{\mathrm{2}h}}\right){\cos }^{\mathrm{2}}\left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right)\right.\\ {\displaystyle }& {\displaystyle }\left..+\mathrm{3}{L}_y\sin \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right)\cos \left({\displaystyle \frac{\mathrm{3}z}{h}}\right){\cos }^{\mathrm{2}}\left({\displaystyle \frac{\mathit{\pi }x}{L_x}}\right)\right)\\ & {\displaystyle }& {\displaystyle }\sin \left({\displaystyle \frac{\mathit{\pi }x}{L_x}}\right)\cos \left({\displaystyle \frac{z}{h}}\right),\end{array}$$ $$\begin{array}{cc}{\displaystyle \frac{\partial \left(wT\right)}{\partial z}}=& {\displaystyle \frac{\mathrm{5}\mathit{\pi }}{L_x{L}_y}}\left[\mathrm{2}{L}_x\left(\cos \left({\displaystyle \frac{z}{\mathrm{2}h}}\right)-\cos \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)\right)\right.\\ {\displaystyle }& {\displaystyle }.\sin \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right)L_y\left(\sin \left({\displaystyle \frac{\mathrm{3}z}{h}}\right)+\sin \left(\mathrm{3}\right)\right)\\ & {\displaystyle }& {\displaystyle }\left..\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{L_x}}x\right)\right]\sin \left({\displaystyle \frac{z}{h}}\right)\sin \left({\displaystyle \frac{\mathit{\pi }x}{L_x}}\right)\sin \left({\displaystyle \frac{\mathit{\pi }y}{L_y}}\right).\end{array}$$ These terms are added as source terms to the right-hand side of Eqs. (), (), (), and (). In the weak form, this corresponds to multiplying the analytical function by the test function and integrating over the domain. The solutions were derived using the SymPy symbolic mathematics Python library .

CPU cost comparison against SLIM

A strong scaling test was carried out with both Thetis and the SLIM 3-D model using the baroclinic eddies test case. These tests were carried out on a Linux cluster with 16-core Intel Xeon E5620 processors and Mellanox Infiniband interconnect. The total time spent to run 40 time steps is presented in Table . The table also lists the speed-up $s_i=T_{\mathrm{0}}/T_i$, where $T_i$ stands for the wall-clock time for $i$ cores, and the parallel efficiency $p=s_i/i$. For an ideal model, the parallel efficiency remains at unity. The results show that on a single core Thetis runs approximately 3.3 times faster than SLIM. On 24 cores, the ratio is 4.0, and on 144 cores, Thetis is still 2.2 times faster.

CPU time in the baroclinic eddies test case for Thetis and the SLIM 3-D model. Both models ran on identical triangular mesh (4 km resolution, 40 vertical levels) using $\mathit{\nu }=\mathrm{20}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ and a 140 s time step. The wall-clock time was recorded over 40 iterations.

No. of	Wall-clock time (s)		Ratio	Speed-up			Efficiency
cores	Thetis	SLIM	$\frac{T_{\mathrm{SLIM}}}{T_{\mathrm{Thetis}}}$	Thetis	SLIM		Thetis	SLIM
1	1778.71	5928.32	3.33	1.00	1.00		1.00	1.00
2	1034.64	4802.34	4.64	1.72	1.23		0.86	0.62
4	500.11	2380.74	4.76	3.56	2.49		0.89	0.62
8	290.61	1284.08	4.42	6.12	4.62		0.77	0.58
16	206.97	675.14	3.26	8.59	8.78		0.54	0.55
20	141.17	524.61	3.72	12.60	11.30		0.63	0.57
24	110.83	440.09	3.97	16.05	13.47		0.67	0.56
32	88.03	330.00	3.75	20.21	17.96		0.63	0.56
40	73.17	260.47	3.56	24.31	22.76		0.61	0.57
48	64.16	222.79	3.47	27.72	26.61		0.58	0.55
64	56.62	158.31	2.80	31.41	37.45		0.49	0.59
80	49.48	127.95	2.59	35.95	46.33		0.45	0.58
96	43.64	109.10	2.50	40.76	54.34		0.42	0.57
112	39.68	95.24	2.40	44.83	62.25		0.40	0.56
128	36.91	83.37	2.26	48.19	71.11		0.38	0.56
144	35.76	78.05	2.18	49.74	75.96		0.35	0.53

TK designed and implemented most of the solver and carried out the numerical simulations. SK and LM contributed to the design and implementation of the model. AB, DH, and MP supervised the work and guided the implementation of the model and the manuscript.

The authors declare that they have no conflict of interest.

Acknowledgements

The National Science Foundation partially supported this research through cooperative agreement OCE-0424602. The National Oceanic and Atmospheric Administration (NA11NOS0120036 and AB-133F-12-SE-2046), Bonneville Power Administration (00062251), and Corps of Engineers (W9127N-12-2-007 and G13PX01212) provided partial motivation and additional support. This work was supported by the UK's Engineering and Physical Science Research Council (grant numbers EP/M011054/1, EP/L000407/1) and the Natural Environment Research Council (grant number NE/K008951/1). This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575. The authors acknowledge the Texas Advanced Computing Center (TACC) at the University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. Edited by: James Annan Reviewed by: James Annan and one anonymous referee