A High-Order Well-Balanced Discontinuous Galerkin Method for Hyperbolic Balance Laws Based on the Gauss-Lobatto Quadrature Rules

Abstract

In this paper, we develop a high-order well-balanced discontinuous Galerkin method for hyperbolic balance laws based on the Gauss-Lobatto quadrature rules. Important applications of the method include preserving the non-hydrostatic equilibria of shallow water equations with non-flat bottom topography and Euler equations in gravitational fields. The well-balanced property is achieved through two essential components. First, the source term is reformulated in a flux-gradient form in the local reference equilibrium state to mimic the true flux gradient in the balance laws. Consequently, the source term integral is discretized using the same approach as the flux integral at Gauss-Lobatto quadrature points, ensuring that the source term is exactly balanced by the flux in equilibrium states. Our method differs from existing well-balanced DG methods for shallow water equations with non-hydrostatic equilibria, particularly in the aspect that it does not require the decomposition of the source term integral. The effectiveness of our method is demonstrated through ample numerical tests.

1 Introduction

The hyperbolic balance law serves as a fundamental tool for characterizing flow and transport phenomena. In one space dimension, it typically takes the form of

$$\begin{aligned} \textbf{u}_t+\textbf{f}(\textbf{u})_x=\textbf{s}(\textbf{u},x), \quad x\in \mathbb {R}, t>0, \end{aligned}$$

(1.1)

where $\textbf{u}\in \mathbb {R}^m$ is the vector of balanced quantities, $\textbf{f}$ the vector of fluxes and $\textbf{s}$ the vector of source terms. Here and henceforth, letters in boldface font are used to denote quantities of vectors.

Due to the presence of the source term, Eq. (1.1) often admits non-trivial steady-state solutions, which are time-independent and exactly balance the flux gradient with the source term, i.e. $\textbf{f}(\textbf{u}_e(x))_x=\textbf{s}(\textbf{u}_e(x),x)$, where $\textbf{u}_e$ is a steady-state solution. It is desirable to design well-balanced numerical methods that satisfy a discrete version of such a balance, as solutions to many real-world problems are minor deviations from steady states and well-balanced schemes are highly efficient in resolving such small deviations.

A prototypical example of hyperbolic balance laws is the system of shallow water equations, which has been widely used in modeling flows in rivers and coastal areas. The shallow water equations with non-flat bottom topography is written as

$$\begin{aligned} {\left\{ \begin{array}{ll} h_t+m_x=0,\\ m_t+\left( h u^2+\frac{1}{2}gh^2\right) _x=-gh\phi _x, \end{array}\right. } \end{aligned}$$

(1.2)

where h is the water depth, $m=hu$ is the discharge, u is the velocity, g is the gravitational constant, and $\phi$ is the prescribed bottom topography.

The shallow water Eq. (1.2) admit the steady state of lake-at-rest, i.e.

$$\begin{aligned} u=0,\quad h+\phi =\text {constant}, \end{aligned}$$

(1.3)

which is hydrostatic as the velocity is zero. The general non-hydrostatic equilibrium of (1.2) is characterized by

$$\begin{aligned} m=\text {constant},\quad \frac{1}{2}u^2+g(h+\phi )=\text {constant}, \end{aligned}$$

(1.4)

where m is not necessarily zero, thereby the lake-at-rest (1.3) is a special case of (1.4).

Numerous numerical methods have been developed for the exact preservation of the state of lake-at-rest (1.3), see, e.g. [2, 3, 6, 28, 41, 52,53,54], and references therein. It is much more challenging yet highly meaningful to design well-balanced schemes for the non-hydrostatic equilibrium (1.4), which naturally covers the steady state (1.3). We refer the readers to works in the non-exhaustive list [4, 5, 9, 10, 26, 36, 39, 40, 42, 51].

Another important example of the hyperbolic balance laws (1.1) is the system of Euler equations for compressible gas dynamics, which takes the form of

$$\begin{aligned} {\left\{ \begin{array}{ll} \rho _t + m_x = 0,\\ m_t + (\rho u^2+p)_x = -\rho \phi _x,\\ E_t + \left( (E+p)u\right) _x = -m\phi _x, \end{array}\right. } \end{aligned}$$

(1.5)

in gravitational fields, where $\rho$ is the density, $m=\rho u$ is the momentum, u is the velocity, $E=\frac{1}{2}\rho u^2+\frac{p}{\gamma -1}$ is the total energy per unit volume, p is the pressure, $\gamma >1$ is the ratio of specific heats, and $\phi$ is the prescribed potential of the gravitational filed.

The complexity of the Euler equations leads to a richer variety of steady states. Specifically, two types of hydrostatic equilibria have been most thoroughly investigated in previous studies, e.g. [8, 20, 24, 27, 33, 47, 49, 55], which are the isothermal hydrostatic equilibrium,

$$\begin{aligned} u=0,\quad \frac{p}{\rho }=\text {constant},\quad \frac{p}{\rho }\log \rho +\phi =\text {constant}, \end{aligned}$$

(1.6)

and the polytropic hydrostatic equilibrium,

$$\begin{aligned} u=0,\quad \frac{p}{\rho ^{\gamma }}=\text {constant},\quad \frac{\gamma }{\gamma -1}\frac{p}{\rho }+\phi =\text {constant}. \end{aligned}$$

(1.7)

The steady state (1.7) is a particular case of the general isentropic equilibria,

$$\begin{aligned} m=\text {constant},\quad \frac{p}{\rho ^{\gamma }}=\text {constant},\quad \frac{1}{2}u^2+\frac{\gamma }{\gamma -1}\frac{p}{\rho }+\phi =\text {constant}, \end{aligned}$$

(1.8)

see the studies in [19, 21]. Due to the break of energy conservation, there is generally no corresponding non-hydrostatic equilibria of (1.6).

The objective of this paper is to establish a high-order, well-balanced discontinuous Galerkin (DG) method for the hyperbolic balance law (1.1), which in particular is capable of preserving the general non-hydrostatic equilibria (1.4) and (1.8) for the shallow water Eq. (1.2) and Euler equations (1.5), respectively. The DG method was first introduced in 1973 by Reed and Hill [46] for solving linear steady-state hyperbolic balance laws, and was later developed into the Runge–Kutta discontinuous Galerkin (RKDG) by Cockburn et al. in a series of papers [11,12,13,14,15] for nonlinear time-dependent problems. The DG method adopts piecewise polynomials to form its function space and incorporates the concept of numerical fluxes from finite volume methods, making it especially well-suited for solving hyperbolic problems due to its advantages in high-order accuracy, compact stencil, flexibility in geometry, and ease of parallelism. There are many researches on well-balanced methods within the DG framework. For instance, Xing and Shu developed a high order well-balanced DG method in [53] for hyperbolic balance laws with separable source terms. Subsequently, a more straightforward well-balanced DG for shallow water equations was proposed in [54] based on the hydrostatic reconstruction [2], thanks to the fact that the hydrostatic equilibrium solution lives in the DG space if the bottom topography is projected onto the same space. However, great difficulties arise when designing well-balanced methods for non-hydrostatic equilibria, as the steady-state solutions generally do not live in the space of piecewise polynomials. In [51], the problem was addressed by a recovery technique of local reference equilibrium states and decomposition of the source term into the equilibrium and residual parts. The method has been extended to preserve the moving-water equilibrium of the Ripa model in [5] and the polytropic equilibrium of Euler equations in [33], and incorporated with a modified oscillation free discontinuous Galerkin (OFDG) damping term in [38] to achieve better robustness. More recently, Zhang et al. [59] proposed a new well-balanced DG scheme for the non-hydrostatic equilibria of shallow water equations. The scheme directly approximates the equilibrium variables in the DG space, rather than the conserved ones, resulting in a much cleaner approach yet with increased computational costs in Newton iteration. We also refer to well-balanced DG methods based on global fluxes in [39, 40], as well as related works cited therein.

In this paper, we propose a novel well-balanced DG method for hyperbolic balance laws based on the Gauss-Lobatto quadrature rules. The method is aimed at preserving the non-hydrostatic equilibria of nodal values at the $k+1$ Gauss-Lobatto points, which can be chosen as the degrees of freedom associated with the Lagrange basis for simplicity, in a $P^k$-DG scheme. The key components of the proposed method are the flux-gradient reformulation of the source terms in local reference equilibrium states and the same approach of discretization for the source term and flux integral at the Gauss-Lobatto points. The special reformulation of the source term is inspired by the treatments in [8, 32, 34], see also [21, 43] for closely related treatment and the survey [23]. A high order well-balanced finite difference WENO scheme has been proposed in our recent work [58] based on the same reformulation. In spirit, our method is very close to the work in [8], which used the nodal DG method to preserve the hydrostatic equilibria of the Euler equations in the gravitational fields.

The remaining of the paper is organized as follows. In Sect. 2, we review some important properties of the steady states of the shallow water equations and Euler equations. Consequently, the well-balanced DG scheme is established in Sect. 3 for these equations. The good performance of the scheme is validated in Sect. 4 by ample numerical tests. We discuss and implement the extension of the algorithm to two dimensions in Sect. 5. Finally, we end up with some concluding remarks in Sect. 6.

2 Steady States of the Shallow Water Equations and Euler Equations

In this section, we review some properties of the non-hydrostatic equilibria of the shallow water equations and Euler equations, and introduce necessary definitions to facilitate the construction of our well-balanced DG method.

2.1 The Shallow Water Equations with Non-flat Bottom Topography

For the shallow water equations (1.2), we introduce the vector of equilibrium variables

$$\begin{aligned} \textbf{v}=\left( m,Q\right) ^T,\quad ~\text {where}~ Q=\frac{1}{2}u^2+g(h+\phi ). \end{aligned}$$

(2.1)

The vector is constant in the non-hydrostatic equilibria (1.4). Moreover, we denote the mapping from the vector of conserved variables $\textbf{u}=(h,m)^T$ to equilibrium ones by $\textbf{v}=V(\mathbf {\textbf{u}},\phi )$. Its inverse mapping is denoted by $\textbf{u}=U(\textbf{v},\phi )$. If $\textbf{u}_e$ is an equilibrium state satisfying (1.4), then $\textbf{u}_e(x)=U(\textbf{v}_e,\phi (x))$ for the constant $\textbf{v}_e=V(\textbf{u}_e(x),\phi (x))$.

For the calculation of $U(\textbf{v},\phi )$, we introduce the Froude number

$$\begin{aligned} \text {Fr}=\frac{|u|}{c}, \quad \text {where}~ c=\sqrt{gh}. \end{aligned}$$

(2.2)

The flow regimes of $\textbf{u}$ are classified as subcritical, critical or supercritical if $\text {Fr}<1,=1$ or $>1$. In each flow regime, the water depth h is uniquely determined by the formula of trigonometric solution of cubic equation [26, 30],

$$\begin{aligned} \begin{aligned} h(m,Q,\phi )=\frac{Q-g\phi }{3g}\left( 1+\cos (\frac{\theta }{3})-\sqrt{3}\sigma \sin (\frac{\theta }{3})\right) , \end{aligned} \end{aligned}$$

(2.3)

where $\theta =\arccos \left( \frac{27}{4}\frac{g^2\,m^2}{\left( Q-g\phi \right) ^3}-1\right)$ and $\sigma ={{\,\textrm{sign}\,}}(\text {Fr-1})$, if the equilibrium variables m, Q are given with $Q-g\phi \ge \frac{3}{2}\left( g|m|\right) ^{\frac{2}{3}}$. For more properties of the equilibrium variables and the mapping $U(\textbf{v},\phi )$, we refer the readers to [42].

2.2 The Euler Equations in Gravitational Fields

For the Euler equations (1.5), we introduce the vector of equilibrium variables

$$\begin{aligned} \textbf{v}=\left( m,s,Q\right) ^T,\quad \text {where}~s=\frac{p}{\rho ^\gamma }, Q=\frac{1}{2}u^2+\frac{\gamma s}{\gamma -1}\rho ^{\gamma -1}+\phi . \end{aligned}$$

(2.4)

The vector is constant in the isentropic non-hydrostatic equilibria (1.8). Similarly, we denote the mapping from the vector of conserved variables $\textbf{u}=(\rho ,m,E)^T$ to equilibrium ones by $\textbf{v}=V(\textbf{u},\phi )$. Its inverse mapping is denoted by $\textbf{u}=U(\textbf{v},\phi )$. If $\textbf{u}_e$ is an equilibrium state satisfying (1.8), then $\textbf{u}_e(x)=U(\textbf{v}_e,\phi (x))$ for the constant $\textbf{v}_e=V(\textbf{u}_e(x),\phi (x))$.

For the calculation of $U(\textbf{v},\phi )$, we introduce the Mach number

$$\begin{aligned} \text {Mach} = \frac{|u|}{c}, \quad \text {where}~c=\sqrt{\gamma p/\rho }. \end{aligned}$$

(2.5)

The flow regimes of $\textbf{u}$ are classified as subcritical, critical or supercritical if $\text {Mach}<1, =1$ or $>1$. In each flow regime, the density $\rho$ is uniquely determined by the equation

$$\begin{aligned} \frac{1}{2}\frac{m^2}{\rho ^2}+\frac{\gamma s}{\gamma -1}\rho ^{\gamma -1}+\phi =Q, \end{aligned}$$

(2.6)

if the equilibrium variables m, s, Q are given with $Q-\phi \ge \left( \frac{1}{2}+\frac{1}{\gamma -1}\right) \left( \gamma s\right) ^{\frac{2}{\gamma +1}}|m|^{\frac{2(\gamma -1)}{\gamma +1}}$. In practice, we can use Newton iteration to calculate $\rho$ from (2.6) and obtain $E=\frac{1}{2}\frac{m^2}{\rho }+\frac{s\rho ^\gamma }{\gamma -1}$ consequently. For more properties of the equilibrium variables and mapping $U(\textbf{v},\phi )$, we refer to [21].

Remark 2.1

The hydrostatic isentropic profiles (i.e., $u=0$ in (1.8) ) in a gravity field are unstable to convection, a physical instability [31, Chapter 1, Section 4]. Convection phenomena occur on the fluid advection time scale, which is usually much slower than acoustic phenomena. Well-balanced schemes are ideally suited to simulate such convection phenomena, as they precisely resolve the perturbations of the equilibrium (i.e., convection) without any spurious numerical errors due to truncation errors stemming from the resolution of the hydrostatic equilibrium. Convection simulations using well-balanced schemes were exemplified, for instance, in [25].

3 The Well-Balanced Discontinuous Galerkin Method

In this section, we establish the well-balanced DG method for the shallow water equations (1.2) and Euler equations (1.5) in a unified framework, and give the proof of its well-balanced property.

3.1 Notations

We first introduce the notations to be used throughout this section. Consider the domain $\Omega =[a,b]\subset \mathbb {R}$, we take the partition $a=x_{\frac{1}{2}}<x_{\frac{3}{2}}<\cdots <x_{N+\frac{1}{2}}=b$, and denote the j-th cell by $I_j=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$ with size $\Delta x_j=x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}$, for $j=1,2,\ldots ,N$.

The function space of the $P^k$-DG scheme is defined as

$$\begin{aligned} W_{h}^{k}=\{w\in L^2([a,b]): w|_{I_j}\in P^k(I_j), j=1,2,\ldots ,N\}, \end{aligned}$$

(3.1)

where $P^k(I)$ denotes the space of polynomials of order no greater than k on the interval I. Since $w\in W_{h}^{k}$ are piecewise polynomials with possible discontinuities at interfaces, we define the left/right limits at $x_{j+\frac{1}{2}}$ by $w_{j+\frac{1}{2}}^{\pm }=\lim _{\epsilon \rightarrow 0^+} w(x_{j+\frac{1}{2}}\pm \epsilon )$. We also denote $\textbf{W}_{h}^{k}=\left[ W_{h}^{k}\right] ^m$ and $\textbf{w}_{j+\frac{1}{2}}^{\pm }=\lim _{\epsilon \rightarrow 0^+} \textbf{w}(x_{j+\frac{1}{2}}\pm \epsilon )$ for $\textbf{w}\in \textbf{W}_{h}^{k}$ as the corresponding vector version of the function space and left/right limits.

To establish the well-balanced DG method, we adopt the $(k+1)$-point Gauss-Lobatto quadrature rule in the $P^k$-DG scheme, and denote the quadrature of a continuous function g on $I_j$ by to distinguish it from the exact integral $\int _{I_j}g dx$, where $\{x_{i}^{j}\}_{i=1}^{k+1}$ and $\{\omega _i\}_{i=1}^{k+1}$ are the Gauss-Lobatto quadrature points (on $I_j$) and weights, respectively. Note that $x_{1}^{j}=x_{j-\frac{1}{2}}$ and $x_{k+1}^{j}=x_{j+\frac{1}{2}}$, which is important for the scheme to achieve the well-balanced property. We denote by $\mathcal {I}^{j}\left[ \cdot \right]$ the Lagrange interpolation operator at the Gauss-Lobatto points on $I_j$, i.e. $\mathcal {I}^j[g]\in P^{k}(I_j)$ with $\mathcal {I}^{j}\left[ g\right] (x_{i}^{j})=g(x_{i}^{j})$ for $i=1,\ldots , k+1$. For convenience, we write the shortened notation $\mathcal {I}[g]$ for $\mathcal {I}^{j}[g]$, if the cell j is clear from the context.

3.2 The Well-Balanced DG Scheme

To construct the scheme, we reformulate the hyperbolic balance law (1.1) on $I_j$ as follows [23],

$$\begin{aligned} \textbf{u}_t+\textbf{f}(\textbf{u})_x=\frac{\textbf{s}(\textbf{u},x)}{\textbf{s}(\textbf{u}_{e}^j,x)}\textbf{f}(\textbf{u}_{e}^j)_x, \end{aligned}$$

(3.2)

where $\textbf{u}_{e}^{j}$ is a local reference equilibrium state satisfying $\textbf{f}(\textbf{u}_{e}^j(x))_x=\textbf{s}(\textbf{u}_{e}^{j}(x),x)$ on $I_j$, whose definition is provided below, and the multiplication and division on the right-hand side must be understood component-wise and set to 1 if a $\frac{0}{0}$ appears in the ratio of the equilibrium flux to the source term.

The well-balanced DG scheme based on the reformulation (3.2) is given as follows: Find $\textbf{U}(t)\in \textbf{W}_{h}^{k}$, such that, $\forall \textbf{w}\in \textbf{W}_{h}^{k}$,

(3.3)

for $j=1,2,\ldots ,N$, where

$$\begin{aligned} \begin{aligned} \hat{\textbf{f}}_{j+\frac{1}{2}}^{l}=\textbf{F}\left( \textbf{U}_{j+\frac{1}{2}}^{*,-},\textbf{U}_{j+\frac{1}{2}}^{*,+}\right) -\textbf{f}(\textbf{U}_{j+\frac{1}{2}}^{*,-})+\textbf{f}(\textbf{U}_{j+\frac{1}{2}}^{-}),\\ \hat{\textbf{f}}_{j-\frac{1}{2}}^{r}=\textbf{F}\left( \textbf{U}_{j-\frac{1}{2}}^{*,-},\textbf{U}_{j-\frac{1}{2}}^{*,+}\right) -\textbf{f}(\textbf{U}_{j-\frac{1}{2}}^{*,+})+\textbf{f}(\textbf{U}_{j-\frac{1}{2}}^{+}), \end{aligned} \end{aligned}$$

(3.4)

are numerical fluxes based on the flux function $\textbf{F}(\cdot ,\cdot )$, e.g. the Lax-Friedrichs flux

$$\begin{aligned} \textbf{F}(\textbf{U}^-, \textbf{U}^+)=\frac{1}{2}\left( \textbf{f}(\textbf{U}^-)+\textbf{f}(\textbf{U}^+)-\alpha \left( \textbf{U}^+-\textbf{U}^-\right) \right) , \end{aligned}$$

(3.5)

with $\alpha =\max _{\textbf{U}}\{|\lambda _1(\textbf{U})|,\ldots ,|\lambda _m(\textbf{U})|\}$ and $\lambda _i$ the eigenvalues of the system (1.1), and the generalized hydrostatic reconstruction [2, 7, 17, 51, 59] $\textbf{U}_{j+\frac{1}{2}}^{*,\pm }$ are defined as

$$\begin{aligned} \textbf{U}_{j+\frac{1}{2}}^{*,\pm }=U(\textbf{V}_{j+\frac{1}{2}}^{\pm },\phi ^{*}_{j+\frac{1}{2}}), \end{aligned}$$

(3.6)

with $\textbf{V}_{j+\frac{1}{2}}^{\pm }=V(\textbf{U}_{j+\frac{1}{2}}^{\pm }, \phi _{j+\frac{1}{2}}^{\pm })$ and $\phi _{j+\frac{1}{2}}^{*}=\min (\phi _{j+\frac{1}{2}}^{-},\phi _{j+\frac{1}{2}}^{+})$. The Lagrange basis at Gauss-Lobatto points is particularly suitable for the computation of the scheme, and the same quadrature can also be used for the left-hand side, which results in the nodal formulation [22].

The local reference equilibrium state $\textbf{u}_e^j$ is defined as follows,

$$\begin{aligned} \textbf{u}_{e}^{j}(x)=U(\textbf{V}_{e}^{j},\phi (x)),\quad x\in I_j, \end{aligned}$$

(3.7)

where $\textbf{V}_{e}^{j}=V(\textbf{U}(x_{\ell }^j),\phi (x_{\ell }^{j}))$, $\ell =\arg \max _{1\le i\le k+1}\phi (x_{i}^{j})$, and the values of piecewise functions at $x_{1}^{j}$ and $x_{k+1}^{j}$ are understood as inside $I_j$ for the time being. With such a choice of $\ell$, $\textbf{u}_{e}^{j}$ is always well-defined on $I_j$.

Remark 3.1

Note that if $\phi _{j+\frac{1}{2}}^{-}=\phi _{j+\frac{1}{2}}^{+}$, we have $\textbf{U}_{j+\frac{1}{2}}^{\star ,\pm }=\textbf{U}_{j+\frac{1}{2}}^{\pm }$, thus the numerical flux $\hat{\textbf{f}}_{j+\frac{1}{2}}^{l}=\hat{\textbf{f}}_{j+\frac{1}{2}}^{r}=\textbf{F}(\textbf{U}_{j+\frac{1}{2}}^{-},\textbf{U}_{j+\frac{1}{2}}^{+})$ is just the standard one when $\phi$ is continuous. The numerical flux in (3.4) is non-standard when $\phi$ has discontinuities. Since the gravitational potential $\phi$ in the Euler equations is typically continuous, as it is a solution of the Poisson equations, we explain the numerical flux (3.4) using the example of the shallow water equations with a discontinuous bottom $\phi$. This explanation also applies to more sophisticated cases, e.g., the shallow water equations in open channels [45] or the Ripa model [16].

We denote by $[\phi ]_{j+\frac{1}{2}}=\phi _{j+\frac{1}{2}}^{+}-\phi _{j+\frac{1}{2}}^{-}$, and define $\tilde{\phi }(s)=(1-s)\phi _{j+\frac{1}{2}}^{-}+s\phi _{j+\frac{1}{2}}^{+},\,s\in [0,1]$, such that $\tilde{\phi }(0)=\phi _{j+\frac{1}{2}}^{-}$, $\tilde{\phi }(1)=\phi _{j+\frac{1}{2}}^{+}$, and $\tilde{\phi }_s(s)=[\phi ]_{j+\frac{1}{2}}$. The following equation holds:

$$\begin{aligned} \textbf{f}(U(\textbf{v},\tilde{\phi }(s)))_s=\textbf{s}(U(\textbf{v},\tilde{\phi }(s))),\quad s\in [0,1] \end{aligned}$$

for fixed $\textbf{v}$, as the steady states of the shallow water equations are characterized by constant equilibrium variables. Without loss of generality, we assume $\phi _{j+\frac{1}{2}}^{-}<\phi _{j+\frac{1}{2}}^{+}$. By definition, we have $\phi _{j+\frac{1}{2}}^{\star }=\phi _{j+\frac{1}{2}}^{-}$, $\textbf{U}_{j+\frac{1}{2}}^{\star ,\pm }=U(\textbf{V}_{j+\frac{1}{2}}^{\pm },\phi _{j+\frac{1}{2}}^-)$, and consequently,

$$\begin{aligned} \begin{aligned} \hat{\textbf{f}}_{j+\frac{1}{2}}^{r}-\hat{\textbf{f}}_{j+\frac{1}{2}}^{l}=&\,\textbf{f}(U(\textbf{V}_{j+\frac{1}{2}}^{+},\phi _{j+\frac{1}{2}}^{+}))-\textbf{f}(U(\textbf{V}_{j+\frac{1}{2}}^{+},\phi _{j+\frac{1}{2}}^{-}))\\ =&\,\textbf{f}(U(\textbf{V}_{j+\frac{1}{2}}^{+},\tilde{\phi }(1)))-\textbf{f}(U(\textbf{V}_{j+\frac{1}{2}}^{+},\tilde{\phi }(0)))\\ =&\int _{0}^{1}\textbf{f}(U(\textbf{V}_{j+\frac{1}{2}}^{+}, \tilde{\phi }(s)))_s \,ds =\int _{0}^{1}\textbf{s}\left( U(\textbf{V}_{j+\frac{1}{2}}^{+},\tilde{\phi }(s))\right) \,ds\\ =&\int _{0}^{1} \begin{pmatrix} 0\\ -gh(\textbf{V}_{j+\frac{1}{2}}^{+},\tilde{\phi }(s))\tilde{\phi }_s(s) \end{pmatrix}ds\\ =&\begin{pmatrix} 0\\ -g\bar{h}_{j+\frac{1}{2}}[\phi ]_{j+\frac{1}{2}} \end{pmatrix}, \end{aligned} \end{aligned}$$

where $\bar{h}_{j+\frac{1}{2}}=\int _{0}^{1} h(\textbf{V}_{j+\frac{1}{2}}^{+},(1-s)\phi _{j+\frac{1}{2}}^{-}+s\phi _{j+\frac{1}{2}}^{+}) \,ds$ is an approximation to $h_{j+\frac{1}{2}}$ in a particular equilibrium path. The first equality is because $\textbf{U}_{j+\frac{1}{2}}^{\pm }=U(\textbf{V}_{j+\frac{1}{2}}^{\pm },\phi _{j+\frac{1}{2}}^{\pm })$ and the function $h(\cdot )$ represents the first component of $U(\cdot )$ in the proof.

Therefore, the difference of the non-conservative fluxes $\hat{\textbf{f}}_{j+\frac{1}{2}}^{r}-\hat{\textbf{f}}_{j+\frac{1}{2}}^{l}$ contributes a Dirac-$\delta$ type source $-g\bar{h}_{j+\frac{1}{2}}[\phi ]_{j+\frac{1}{2}}\delta (x-x_{j+\frac{1}{2}})$ at the interface $x_{j+\frac{1}{2}}$ to the momentum equation $m_t+(hu^2+\frac{1}{2}gh^2)_x=-gh\phi _{x}$. We refer the readers to the original papers [7, 17] for further explanations of the fluxes (3.4).

Remark 3.2

For the Euler equations, we have $\textbf{s}(\textbf{u},x)=(0,-\rho \phi _x,-m\phi _x)^T, \textbf{s}(\textbf{u}_{e}^{j},x)=(0,-\rho _{e}^{j}\phi _x,-m_{e}^{j}\phi _x)^T$, and $\textbf{f}(\textbf{u}_{e}^{j})=(m_{e}^{j},\frac{(m_{e}^{j})^2}{\rho _{e}^{j}}+p_{e}^{j}, (E_{e}^{j}+p_{e}^{j})\frac{m_{e}^{j}}{\rho _{e}^{j}})$, where $m_{e}^{j}$ is a constant. After cancellations, the reformulation of the source term in (3.2) can be simplified as $\left( 0, \frac{\rho }{\rho _{e}^{j}}, m\right) ^T\times \partial _x(0,\frac{(m_{e}^{j})^2}{\rho _{e}^{j}}+p_{e}^{j}, (E_{e}^{j}+p_{e}^{j})\frac{1}{\rho _{e}^{j}})^T$, where the operations are understood component-wisely. Similarly, for the shallow water equations, the reformulation of the source term in (3.2) can be simplified as $\left( 0,\frac{h}{h_{e}^{j}}\right) ^T\times \partial _x(0,\frac{(m_{e}^{j})^2}{h_{e}^{j}}+\frac{1}{2}g(h_{e}^{j})^2)^T$.

3.3 The Well-Balanced Property

The scheme (3.3) is well-balanced in the sense that, if its solution $\textbf{U}$ is in equilibrium at all Gauss-Lobatto points $\{x_{i}^{j}\}_{i=1}^{k+1}, j=1,2,\ldots ,N$, then the right-hand side is zero, which can be stated as the following theorem.

Theorem 3.1

If $\textbf{U}\in \textbf{W}_{h}^{k}$ satisfies $\textbf{U}|_{I_j}=\mathcal {I}^j\left[ \textbf{u}_e\right] , j=1,\ldots ,N$, for an equilibrium state $\textbf{u}_{e}$ of (1.1), then the right-hand side of (3.3) is zero.

Proof

We denote by $\textbf{V}_{i}^{j}=V(\textbf{U}(x_{i}^{j}),\phi (x_{i}^{j}))$ for $i=1,\ldots ,k+1, j=1,\ldots ,N$. Since $\textbf{U}(x_{i}^{j})=\textbf{u}_e(x_{i}^{j})$ and $V(\textbf{u}_e(x),\phi (x))\equiv \textbf{v}_e$, we have $\textbf{V}_{i}^{j}\equiv \textbf{v}_e$, $\forall i,j$, where $\textbf{v}_e$ is a constant. Therefore, $\textbf{U}_{j+\frac{1}{2}}^{*,-}=U(\textbf{v}_e,\phi _{j+\frac{1}{2}}^{*})=\textbf{U}_{j+\frac{1}{2}}^{*,+}$, which implies $\hat{\textbf{f}}_{j+\frac{1}{2}}^{l}=\textbf{f}(\textbf{U}_{j+\frac{1}{2}}^{-})$ and $\hat{\textbf{f}}_{j-\frac{1}{2}}^{r}=\textbf{f}(\textbf{U}_{j-\frac{1}{2}}^{+})$ from the consistency of $\textbf{F}(\cdot ,\cdot )$. Moreover, by the definition (3.7), we have $\textbf{u}_{e}^{j}(x)=\textbf{u}_e, j=1,\ldots , N$.

Consequently, the right-hand side of (3.3) is

where the second equality is because the $(k+1)$-point Gauss-Lobatto quadrature rule is exact for polynomials of degree no greater than $2k-1$, and the last equality is a consequence of integration by parts. $\square$

Remark 3.3

By Theorem 3.1, if the initial condition $\textbf{u}_0$ of Eq. (1.1) is in equilibrium, which is generally non-polynomial, we should take the initial condition $\textbf{U}_0$ of scheme (3.3) by the Lagrange interpolation $\textbf{U}_0=\mathcal {I}\left[ \textbf{u}_0\right]$ to preserve the equilibrium exactly.

3.4 An Extra Damping Term

Following [38], we add an extra oscillation free discontinuous Galerkin (OFDG) damping term to the scheme (3.3) to control the spurious oscillations without breaking the well-balanced property.

The well-balanced OFDG scheme is given as follows: Find $\textbf{U}(t)\in \textbf{W}_{h}^{k}$, such that, $\forall \textbf{w}\in \textbf{W}_{h}^{k}$,

$$\begin{aligned} LHS_j=RHS_j-\sum _{\ell =0}^{k}\frac{\sigma _{j}^{\ell }(\textbf{U})}{\Delta x_j}\int _{I_j}\left( \textbf{U}_{r}^{j}-P_{h}^{\ell -1}\textbf{U}_{r}^{j}\right) \cdot \textbf{w}dx, \end{aligned}$$

(3.8)

for $j=1,2,\ldots ,N$, where $LHS_j$ and $RHS_j$ are the original left-hand and right-hand sides of (3.3), respectively, and $\textbf{U}_{r}^{j}=\textbf{U}|_{I_j}-\mathcal {I}\left[ \textbf{u}_{e}^{j}\right]$ is the residual of equilibrium state, $P_{h}^{\ell -1}, \ell \ge 1$ is the standard component-wise $L^2$-projection on $P^{\ell -1}(I_j)$ and $P^{-1}_{h}=P_{h}^{0}$. The damping coefficients $\sigma _{j}^{\ell }(\textbf{U})$ are constructed in exactly the same way as [38] and [37] for the shallow water equations and Euler equations, respectively. Since $\textbf{U}|_{I_j}=\mathcal {I}[\textbf{u}_{e}^{j}]$ in equilibrium, it is clear that the OFDG damping term does not destroy the well-balanced property of (3.3).

Beside the OFDG method, appropriate slope limiters, e.g., [5, 59], may also be used to control the oscillation without affecting the well-balanced property.

4 Numerical Tests

In this section, we test the performance of the proposed well-balanced DG scheme (3.8) for the shallow water equations and Euler equations, with the classic fourth order Runge-Kutta time discretization. We take the step size of time as $\Delta t=\frac{\text {CFL}}{\left| \left| |u|+c\right| \right| _{L^\infty (\Omega )}}\Delta x$, where $\Delta x$ is the minimum cell length of mesh and CFL is a constant taken as 0.1 in the tests if not otherwise stated. The gravitational constant in the shallow water equations are taken as $g=9.812$ in all tests. To save space, we only show the results of the $P^2$-DG scheme in the well-balanced tests. All numerical examples in this section are classic. We refer the readers to the results of other methods, e.g. [10, 21, 29, 36, 42, 48, 55, 59], for comparison.

4.1 The Shallow Water Equations

In this subsection, we test the well-balanced scheme for the shallow water equations.

Example 4.1

Accuracy test In this example, we test the accuracy of the proposed well-balanced DG method for the shallow water equations with a smooth bottom.

The bottom topography and initial conditions of the problem are given by

$$\begin{aligned} \phi (x)=\sin ^2(\pi x),\quad h(x,0)=5+e^{\cos (2\pi x)},\quad m(x,0)=\sin (\cos (2\pi x)). \end{aligned}$$

(4.1)

The computational domain is $\Omega =[0,1]$ with periodic boundary condition.

Since the analytical solution of the problem is not available, we calculate the reference solution using the ninth-order WENO finite difference scheme on $N=6400$ grid. The $L^1$ errors and orders of convergence of h and m with N uniform cells are listed in Table 1, at the terminal time $T=0.1$ before shock formation. We can clearly observe high order convergence rates for $P^1$, $P^2$ and $P^3$-DG methods from the table.

Table 1 $L^1$ errors and orders of convergence in the accuracy test 4.1 for the shallow water equations

Example 4.2

Well-balanced test for a continuous bottom topography In this example, we test the well-balanced property of our scheme for a continuous bottom in different flow regimes.

The bottom topography of the problem is given by

$$\begin{aligned} \phi (x)= {\left\{ \begin{array}{ll} 0.2-0.05(x-10)^2,\quad & 8\le x\le 12,\\ 0,& \text {otherwise}, \end{array}\right. } \end{aligned}$$

(4.2)

on the domain $\Omega =[0,25]$. We impose the boundary conditions by keeping the values on boundaries constant and equal to the initial conditions.

Below, we consider four non-hydrostatic equilibria in different flow regimes:

(a) Subcritical flow:

$$\begin{aligned} m(x,0)=4.42,\quad Q(x,0)=22.06605, \quad \text {with}\quad \text {Fr}(x,0)<1. \end{aligned}$$

(b) Supercritical flow:

$$\begin{aligned} m(x,0)=24,\quad Q(x,0)=91.624, \quad \text {with}\quad \text {Fr}(x,0)>1. \end{aligned}$$

(c) Transcritical flow without shock:

$$\begin{aligned} \begin{aligned} m(x,0)=1.53,&\quad Q(x,0)=\frac{3}{2}\left( 9.812\times 1.53\right) ^{\frac{2}{3}}+9.812\times 0.2,\\ \text {with}&\quad {\left\{ \begin{array}{ll} \text {Fr}(x,0)<1, & x<10,\\ \text {Fr}(x,0)>1, & x>10. \end{array}\right. } \end{aligned} \end{aligned}$$

(d) Transcritical flow with a stationary shock:

$$\begin{aligned} \begin{aligned}&m(x,0)=0.18,\quad \\&Q(x,0)= {\left\{ \begin{array}{ll} \frac{3}{2}\left( 9.812\times 0.18\right) ^{\frac{2}{3}}+9.812\times 0.2,\quad & x\le 11.665504281554291,\\ \frac{0.18^2}{2\times 0.33^2}+9.812\times 0.33, & \text {otherwise}, \end{array}\right. } \\&\quad \text {with}\quad {\left\{ \begin{array}{ll} \text {Fr}(x,0)<1, & x<10 ~\text {or}~x>11.665504281554291,\\ \text {Fr}(x,0)>1, & 10<x<11.665504281554291.\\ \end{array}\right. } \end{aligned} \end{aligned}$$

We compute the solutions of case (a), (b), (c) and (d) on $N=200$ meshes up to $T=10$. The cells are uniform for the first three cases. Due to the appearance of a shock at $x=11.665504281554291$ in the last case, we use the Roe flux as in [59] and partition 100 uniform cells for the regions on the left and right sides of the shock, respectively, such that the shock locates exactly at a cell interface of the mesh.

The $L^1$ and $L^\infty$ errors of the numerical solutions compared with the solution at initial time are listed in Table 2, from which we can observe that the errors are in the round-off level, which confirms the well-balanced property of the scheme.

Table 2 $L^1$ and $L^\infty$ errors for different cases in the well-balanced test 4.2 with $N=200$ at $T=10$

Example 4.3

Well-balanced test for a discontinuous bottom topography In this example, we test the well-balanced property of the proposed scheme for a discontinuous bottom in different flow regimes.

The bottom topography of the problem is given by

$$\begin{aligned} \phi (x)= {\left\{ \begin{array}{ll} 0.2,\quad & 8\le x\le 12,\\ 0,& \text {otherwise}, \end{array}\right. } \end{aligned}$$

(4.3)

on the domain $\Omega =[0,25]$. We impose the boundary conditions by keeping the values on boundaries constant and equal to the initial conditions.

Below, we consider three non-hydrostatic equilibria in different flow regimes:

(a) Subcritical flow: The settings are the same as case (a) in Example 4.2.

(b) Supercritical flow: The settings are the same as case (b) in Example 4.2.

(c) Transcritical flow:

$$\begin{aligned} \begin{aligned} m(x,0)=1.53,&\quad Q(x,0)=\frac{3}{2}\left( 9.812\times 1.53\right) ^{\frac{2}{3}}+9.812\times 0.2,\\ \text {with}&\quad {\left\{ \begin{array}{ll} \text {Fr}(x,0)<1, & x<8,\\ \text {Fr}(x,0)=1, & 8 \le x\le 12,\\ \text {Fr}(x,0)>1, & x>12. \end{array}\right. } \end{aligned} \end{aligned}$$

We compute the solutions of case (a), (b) and (c) on $N=200$ uniform meshes up to $T=10$.

The $L^1$ and $L^\infty$ errors of the numerical solutions compared with the solution at initial time are listed in Table 3, from which we can observe that the errors are in the round-off level, which verifies the well-balanced property of the scheme.

Table 3 $L^1$ and $L^\infty$ errors for different cases in the well-balanced test 4.3 with $N=200$ at $T=10$

Example 4.4

Small perturbation of steady states for the continuous bottom In this example, we test the capability of the proposed scheme to capture small perturbations of equilibria for continuous bottom topography. We add $\delta =0.05$ in the region $x\in [5.75, 6.25]$ on top of the initial water depth of case (a), (b), (c) and (d) in Example 4.2, and compute the solutions up to $T=1.5, 1.0, 1.5$ and 3.0, respectively.

We show the numerical results with locally zoomed boxes in Fig. 1, from which we can see the small perturbations are captured very well. The high-resolution reference solutions obtained from an $N=1000$ grid appear to show some small over- and under-shoots in the figures. They can be effectively eliminated by choosing larger damping coefficients (e.g., using coefficients 10 times larger, according to our tests) in the OFDG method. There are possibly two reasons why the current coefficients produce some small over- and under-shoots. First, the scale of the perturbations in the tests is too small to be stabilized by the normal coefficients proposed in the paper [37, 38]. Second, since the damping terms are modified for the purpose of well-balancing, the original damping coefficients are not the optimal ones. However, since the main focus of this paper is not on stabilization techniques, we do not further explore the optimal coefficients for OFDG damping terms and continue to use the original ones proposed in [37, 38].

Fig. 1

Surface level $h+\phi$ and bottom in the test of small perturbations of steady states for continuous bottom in Example 4.4

Example 4.5

Small perturbation of steady states for the discontinuous bottom In this example, we test the capability of the proposed scheme to capture small perturbations of equilibria for discontinuous bottom topography. We add $\delta =0.05$ in the region $x\in [5.75, 6.25]$ on top of the initial water depth of case (a), (b) and (c) in Example 4.3, and compute the solutions up to $T=1.5, 1.0$, and 1.5, respectively.

We show the numerical results with locally zoomed boxes in Fig. 2, from which we can see that the small perturbations are captured very well.

Fig. 2

Surface level $h+\phi$ and bottom in the test of small perturbations of steady states for discontinuous bottom in Example 4.5

Example 4.6

Dam breaking over a rectangular bump In this example, we test the dam breaking over a rectangular bump.

The bottom topography and initial conditions of the problem are given by

$$\begin{aligned} b(x)={\left\{ \begin{array}{ll} 8, & |x-750|\le 187.5,\\ 0,& \text {otherwise}, \end{array}\right. } \quad h(x,0)={\left\{ \begin{array}{ll} 20-b(x),& x\le 750,\\ 15-b(x), & \text {otherwise}, \end{array}\right. }\quad m(x,0)=0. \end{aligned}$$

The computational domain is $\Omega =[0,1500]$. We impose the boundary conditions by keeping the boundary values constant and equal to the initial conditions.

We compute the solution of the DG scheme on $N=200, 400$ uniform meshes with $\text {CFL}=0.05$, and draw the cell averages of surface level $h+\phi$ at $T=15, 60$ in Fig. 3. The reference solutions in the figure are given by WENO finite difference scheme on $N=6400$ grid. The results exhibit good resolution and essentially oscillation free fashion as in [38].

Fig. 3

Surface level $h+\phi$ of the dam breaking problem 4.6

4.2 The Euler Equations

In this subsection, we test the well-balanced scheme for the Euler equations. The transformation between the conserved variables and equilibrium variables in the scheme requires Newton method in root finding. We set the threshold for the stopping criteria of the Newton method at $1000\epsilon$, where $\epsilon$ represents the machine precision, e.g., $\epsilon \approx 2.2\times 10^{-16}$ for double precision. According to our experience, the algorithm is not sensitive to the stopping criteria, and the number of iterations is less than 10 iterations on average. Since we do not need the generalized hydrostatic reconstruction due to the continuity of $\phi$, we perform root finding $(k+1)$ times (at each Gauss-Lobatto point) per cell. Compared with identifying the local reference equilibrium from cell averages [42, 51], recovering the equilibrium from point values offers significant computational advantages. The percentage of computational time is taken less than $15\%$ in all tests (on the Julia platform). We would like to note that there is indeed no need to calculate the reference steady state at each time stage, although we did so in our implementation. We can reduce the computational cost of root finding by performing the calculation only once per time step, which includes four stages in the fourth-order Runge–Kutta methods, or even performing the calculation per several time time steps.

Example 4.7

Accuracy test In this example, we test the accuracy of the proposed well-balanced DG method for the Euler equations in a gravitational field.

We take the linear gravitational field $\phi (x)=x$ and ratio of specific heat $\gamma =\frac{5}{3}$. An exact solution of the Euler equations (1.5) is then given by

$$\begin{aligned} & \rho (x,t)=1+\frac{1}{5}\sin (\pi (x-u_0 t)),\quad u(x,t)=u_0,\quad \nonumber \\ & p(x,t)=p_0+u_0 t-x+\frac{1}{5\pi }\cos (\pi (x-u_0 t)), \end{aligned}$$

(4.4)

in the domain $\Omega =[0,2]$, where $u_0$, $p_0$ are constants taken as $u_0=1, p_0=4.5$ in the computation.

The $L^1$ errors and orders of convergence of $\rho , m$ and E are listed in Table 4 at the terminal time $T=0.1$, where we can clearly observe the optimal convergence rates for $P^1$, $P^2$ and $P^3$-DG methods.

Table 4 $L^1$ errors and orders of convergence in the accuracy test 4.7 for the Euler equations

Example 4.8

Well-balanced test In this example, we test the well-balanced property of the proposed scheme for steady-state isentropic flows.

We take the gravitational filed $\phi (x)=x$ and ratio of specific heat $\gamma =\frac{5}{3}$, with the initial conditions satisfying

$$\begin{aligned} s(x,0)=1,\quad m(x,0)=-M\gamma ^{\frac{1}{2}}, \quad Q(x,0)=\frac{1}{2} M^2\gamma +\frac{\gamma }{\gamma -1}, \end{aligned}$$

(4.5)

on the domain $\Omega =[0,2]$. The boundary conditions are imposed by keeping constant values that equal the initial conditions.

Below, we consider three equilibria in different flow regimes:

(a) Hydrostatic flow: $M=0$.

(b) Subcritical flow: $M=0.01$.

(c) Supercritical flow: $M=2.5$.

We compute the solutions up to $t=10$ on $N=200$ uniform meshes, and gather the $L^1$ and $L^\infty$ errors compared with the initial conditions in Table 5, from which we can observe the steady states are exactly preserved up to round-off errors.

Table 5 $L^1$ and $L^\infty$ errors for different cases in the well-balanced test 4.8 with $N=200$ at $T=10$

Example 4.9

Small perturbation of steady states In this example, we test the capability of the scheme to capture small perturbations of equilibria.

We add $\delta p=A\exp (-100(x-\bar{x})^2)$ on top of the initial pressure of case (a), (b) and (c) in Example 4.8, and keep density and velocity the same as before, where $A=10^{-6}$ and

$$\begin{aligned} \bar{x}={\left\{ \begin{array}{ll} 1.0,& \text {case (a)},\\ 1.1, & \text {case (b)},\\ 1.5, & \text {case (c)}. \end{array}\right. } \end{aligned}$$

We compute the solutions up to $t=0.45, 0.45$ and 0.25 for the case (a), (b) and (c), respectively. The discrepancy of pressure and velocity with respect to their corresponding base steady states are shown in Fig. 4, together with the results of a non-well-balanced DG scheme with the straightforward approximation of source term and normal OFDG damping term for comparison. We can see the well-balanced scheme captures small perturbations much better than the non-well-balanced scheme on coarse meshes.

Fig. 4

Discrepancy of pressure and velocity with respect to their base steady states in small perturbation test 4.9. WB, well-balanced scheme; NWB, non-well-balanced scheme

Example 4.10

Discontinuous wave propogation In this example, we test the capability of the proposed scheme to capture shocks and large gradients of solutions. We enlarge the perturbation in Example 4.9 by taking $A=1$ in the pressure deviation $\delta p$ therein, with all other parameters kept the same.

The numerical results of the well-balanced scheme are shown in Fig. 5, and are compared with those of the non-well-balanced scheme with a straightforward approximation of source terms and normal OFDG damping term. The reference solutions are computed by the WENO finite difference scheme on $N=6400$ mesh. We can see the well-balanced scheme and non-well-balanced scheme produce indistinguishable solutions with essentially oscillation free fashion as expected, and agree very well with those of the WENO schemes on the fine mesh. We noticed two times larger OFDG damping coefficients than the ones proposed in [37] is needed in the computation of case (b), for both the well-balanced and non-well-balanced DG schemes, to avoid negative pressure. This is possibly because of the low pressure near the right boundary.

Fig. 5

Pressure and velocity in discontinuous wave propogation test 4.10

5 Extension to Two Dimensions

Although the primary focus of this paper is on one spatial dimension, the algorithm can be extended to higher dimensions by following the same principles. The extension will be discussed in this section.

We consider the two-dimensional hyperbolic balance laws of the form:

$$\begin{aligned} \textbf{u}_{t}+\textbf{f}(\textbf{u})_x+\textbf{g}(\textbf{u})_{y}=\textbf{s}(\textbf{u},x,y),\quad (x,y)\in \mathbb {R}^2, t>0, \end{aligned}$$

(5.1)

where $\textbf{u}\in \mathbb {R}^{m}$ represents the vector of balanced quantities, $\textbf{f}$ and $\textbf{g}$ are the vectors of fluxes in the x and y directions, respectively, and $\textbf{s}$ is the vector of source terms. In particular, we primarily focus on the shallow water equations with non-flat bottom topography in two spatial dimensions as a prototypical example:

$$\begin{aligned} {\left\{ \begin{array}{ll} h_t+m_x+n_y=0,\\ m_t+(hu^2+\frac{1}{2} gh^2)_x+(huv)_y=-gh\phi _x,\\ n_t+(huv)_x+(hv^2+\frac{1}{2} gh^2)_y=-gh\phi _y,\\ \end{array}\right. } \end{aligned}$$

(5.2)

where h is the water depth, $m=hu$ and $n=hv$ are the discharges in the x and y directions, respectively, (u, v) is the velocity field of the fluid, g is the gravitational constant, and $\phi$ is the prescribed bottom topography.

There are very few studies in the literature on the preservation of genuinely two-dimensional moving-water equilibria. In [42], a pseudo-2D case, where all variables are constant ($v=0$) along the y-direction, is considered on rectangular grids. One of the main difficulties in designing well-balanced methods for the two-dimensional shallow water equations, and more generally for hyperbolic balance laws, is achieving the target steady state. Unlike the one-dimensional case, where the family of equilibrium states is typically described by simple algebraic equations as demonstrated in Sect. 2, the moving-water equilibria in two dimensions require solving partial differential equations, which allow for a multitude of steady states [42]. As a consequence, identifying the reference equilibrium state $\textbf{u}_{e}$ needed in the reformulation of source terms for general non-hydrostatic scenarios is a highly non-trivial task and, in most scenarios, requires numerical computation [1, 35, 50, 56, 57]. For simplicity, we assume the target steady state a priori in the subsequent discussion, and leave the inference of the steady state at run-time as future work. Another difficulty arises from applying generalized hydrostatic reconstruction (GHR) in two dimensions. In one-dimensional contexts, GHR [7, 17] is well-established for addressing the discontinuities of bottoms at cell interfaces and for preserving the moving-water equilibria in shallow water equations. However, in two-dimensional scenarios, research on GHR remains underdeveloped. To circumvent the complexities of GHR in two dimensions, we treat discontinuous bottoms as single-valued functions and interpolate them continuously at Gauss-Lobatto type nodes. Consequently, the a priori steady states are also single-valued and continuously interpolated. This straightforward approach has proven to be well-balanced and stable in our numerical tests.

The remainder of this section is devoted to constructing and proving the well-balanced DG method, as well as its numerical validation.

5.1 The Well-Balanced DG Method

We first introduce the notations to be used throughout this subsection. Consider a two-dimensional bounded domain $\Omega \subset \mathbb {R}^2$ with a triangulation $\mathcal {T}_h=\{K\}$. The function space of the $P^k$-DG method is defined as

$$\begin{aligned} W_{h}^{k}=\{w\in L^2(\Omega ): w|_{K}\in P^k(K), \,\forall K\in \mathcal {T}_h\}, \end{aligned}$$

(5.3)

where $P^k(K)$ represents the space of polynomials of degree no greater than k on the cell K. As before, we use $\textbf{W}_{h}^{k}=[W_{h}^{k}]^m$ for the corresponding vector version of the function space.

To establish the well-balanced DG method, we introduce the Gauss-Lobatto type nodes $\{(x_{i}^{K},y_{i}^{K})\}_{i=1}^{N_k}$ on each triangle K, with $N_k=\frac{(k+1)(k+2)}{2}$ representing the degrees of freedom of the $P^k$ element. These nodes coincide with the $(k+1)$ Gauss-Lobatto points on each edge. Figure 6 illustrates their distribution on equilateral triangles, and their distribution on general triangles follows the affine mapping. The barycentric coordinates of these nodes for $k=1,\ldots ,5$ are provided in Appendix A. For details of the construction and computation of these nodes, we refer the reader to the monograph [22, Chapter 6, Section 1] of nodal DG methods.

Fig. 6

Distribution of the Gauss-Lobatto type nodes (blue dots) on an equilateral triangle K. These nodes coincide with the $(k+1)$ Gauss-Lobatto points on each edges

We then introduce the Lagrange interpolation operator $\mathcal {I}^{K}$ at the Gauss-Lobatto type nodes, defined such that $\mathcal {I}^{K}[g]\in P^{k}(K)$ and $\mathcal {I}^K [g](x_{i}^{K},y_{i}^{K})=g(x_{i}^{K},y_{i}^{K})$ for $i=1,2,\ldots , N_k$, where $g\in C^0(K)$. Moreover, we define the one-dimensional Lagrange interpolation operator $\mathcal {I}^{\partial K}$ at the Gauss-Lobatto points on each edge of $\partial K$ for functions $g\in C^0(\partial K)$. For convenience, we use the abbreviated notation $\mathcal {I}[g]$ to represent $\mathcal {I}^K[g]$ and $\mathcal {I}^{\partial K}[g]$, when the context clearly indicates its meaning.

The well-balanced DG scheme is defined as follows: Find $\textbf{U}(t)\in \textbf{W}_{h}^{k}$ such that for all $\textbf{w}\in \textbf{W}_{h}^{k}$,

$$\begin{aligned} \begin{aligned} \int _{K}\textbf{U}_t\cdot \textbf{w}\,dxdy&=\int _{K}\left( \mathcal {I}\left[ \textbf{f}(\textbf{U})\right] \cdot \textbf{w}_x+\mathcal {I}\left[ \textbf{g}(\textbf{U})\right] \cdot \textbf{w}_y\right) \,dxdy\\ &\quad -\int _{\partial K}\mathcal {I}\left[ {\textbf{H}}(\textbf{U}^{\text {int}},\textbf{U}^{\text {ext}},\textbf{n})\right] \cdot \textbf{w}\, ds\\&\quad +\int _{K}\mathcal {I}\left[ \frac{\textbf{s}(\textbf{U},x,y)}{\textbf{s}(\textbf{u}_e,x,y)}\left( \mathcal {I}[\textbf{f}(\textbf{u}_e)]_x+\mathcal {I}[\textbf{g}(\textbf{u}_e)]_y\right) \right] \cdot \textbf{w}\,dxdy\\ \end{aligned} \end{aligned}$$

(5.4)

where $\textbf{H}$ is a consistent numerical flux function, e.g., the Lax-Friedrichs flux

$$\begin{aligned} \begin{aligned} \textbf{H}(\textbf{U}^{\text {int}},\textbf{U}^{\text {ext}},\textbf{n})&=\frac{1}{2}\left( \textbf{f}(\textbf{U}^{\text {int}})\textbf{n}_1+\textbf{g}(\textbf{U}^{\text {int}})\textbf{n}_2 + \textbf{f}(\textbf{U}^{\text {ext}})\textbf{n}_1+\textbf{g}(\textbf{U}^{\text {ext}})\textbf{n}_2\right) \\ &\quad -\frac{\alpha }{2}\left( \textbf{U}^{\text {ext}} - \textbf{U}^{\text {int}} \right) , \end{aligned} \end{aligned}$$

(5.5)

with $\alpha = \max _{1\le i\le m,\forall \textbf{U}}|\lambda _i(\textbf{U})|$, where $\lambda _i(\textbf{U})$ the i-th eigenvalue of the Jacobian matrix $\textbf{f}'(\textbf{U})\textbf{n}_1+\textbf{g}'(\textbf{U})\textbf{n}_2$, and $\textbf{n}=(\textbf{n}_1,\textbf{n}_2)^T$ is the unit outer normal on $\partial K$. $\textbf{U}^{\text {int}}$ and $\textbf{U}^{\text {ext}}$ are the numerical solutions on the interior and exterior sides of K, respectively. $\textbf{u}_e$ is a known a priori steady state of (5.1), i.e., $\textbf{f}(\textbf{u}_e)_x+\textbf{g}(\textbf{u}_e)_{y}=\textbf{s}(\textbf{u}_e,x,y)$. Similar to Sect. 3, the OFDG damping [37, 38], slope limiters [5, 59], or filters [44], can be adapted to control spurious oscillations without affecting the well-balanced property.

We have the following important result for the scheme (5.4): the scheme is well-balanced in the sense that if its solution $\textbf{U}$ is in equilibrium at all Gauss-Lobatto type points $\{(x_{i}^{K}, y_{i}^{K})\}_{i=1}^{N_k}$ for $K\in \mathcal {T}_h$, then the right-hand side of the scheme is zero. This is stated in the following theorem.

Theorem 5.1

If $\textbf{U}\in \textbf{W}_{h}^{k}$ satisfies $\textbf{U}|_{K}=\mathcal {I}^K\left[ \textbf{u}_e\right] , \forall K\in \mathcal {T}_h$, for an equilibrium state $\textbf{u}_{e}$ of (5.1), then the right-hand side of (5.4) is zero.

Proof

Due to the consistency of the flux function $\textbf{H}$ and $\textbf{U}|_{K}=\mathcal {I}^K\left[ \textbf{u}_e\right]$, we have

$$\begin{aligned} {\textbf{H}}(\textbf{U}^{\text {int}},\textbf{U}^{\text {ext}},\textbf{n})=\textbf{f}(\textbf{u}_e)\textbf{n}_1+\textbf{g}(\textbf{u}_e)\textbf{n}_2, \end{aligned}$$

at the Gauss-Lobatto points on each edge.

Consequently, the right-hand side of (5.4) is

$$\begin{aligned} \begin{aligned} RHS_K=&\int _{K}\left( \mathcal {I}[\textbf{f}(\textbf{u}_e)]\cdot \textbf{w}_x + \mathcal {I}[\textbf{g}(\textbf{u}_e)]\cdot \textbf{w}_y \right) \,dxdy-\int _{\partial K}\mathcal {I}\left[ \textbf{f}(\textbf{u}_e)\textbf{n}_1+\textbf{g}(\textbf{u}_e)\textbf{n}_2\right] \cdot \textbf{w}\,ds\\&+\int _{K}\left( \mathcal {I}[\textbf{f}(\textbf{u}_e)]_x+\mathcal {I}[\textbf{g}(\textbf{u}_e)]_y\right) \cdot \textbf{w}\,dxdy\\ =&0, \end{aligned} \end{aligned}$$

where the first equality is due to $\textbf{U}=\textbf{u}_{e}$ at all Gauss-Lobatto type nodes, and the second equality results from integration by parts. $\square$

Remark 5.1

It is worth commenting on the conservation property of the two-dimensional well-balanced DG scheme (5.4). Unlike the one-dimensional case, where the discharge $m_e$ at equilibrium is constant, in two spatial dimensions, the discharges $(m_e,n_e)^T$ at equilibrium are generally not in the DG function space because divergence-free vector fields are not necessarily polynomials. As a consequence, the first component of the source term in the two-dimensional scheme (5.4) is not identically zero.

However, this does not compromise the conservation property of the scheme. To illustrate, consider the first equation of the scheme (5.4) for the shallow water equations:

$$\begin{aligned} \begin{aligned} \int _{K}h_t w\,dxdy&=\int _{K}\left( m w_x + n w_y \right) \,dxdy - \int _{\partial K}\mathcal {I}[\hat{h}]w\,ds\\&\quad +\int _{K} \left( \mathcal {I}[m_e]_x+ \mathcal {I}[n_e]_y\right) w \,dxdy, \end{aligned} \end{aligned}$$

where $\hat{h}$ is the Lax-Friedrichs flux for water depth.

Mass conservation becomes evident when we take $w=1$, yielding:

$$\begin{aligned} \begin{aligned} \begin{aligned} \int _{K}h_t\,dxdy=&- \int _{\partial K}\mathcal {I}[\hat{h}]\,ds+\int _{K} \left( \mathcal {I}[m_e]_x + \mathcal {I}[n_e]_y\right) \,dxdy\\ =&-\int _{\partial K}\mathcal {I}[\check{h}]\,ds, \end{aligned} \end{aligned} \end{aligned}$$

(5.6)

where $\check{h}=\hat{h}-m_e\textbf{n}_1-n_e\textbf{n}_2$ is a conservative (i.e., single-valued) numerical flux, with $\int _{\partial K}\hat{h}ds = \int _{\partial {K}}\check{h}ds$.

Remark 5.2

The still-water equilibrium (lake-at-rest):

$$\begin{aligned} h+\phi =\text {constant},\quad u=v=0, \end{aligned}$$

(5.7)

is the most basic steady state for the two-dimensional shallow water equations. For such a particular case, the following local reference equilibrium is available on each cell K:

$$\begin{aligned} \textbf{u}_{e}^{K}= \begin{pmatrix} h_{e}^{K}\\ 0\\ 0 \end{pmatrix},\quad h_{e}^{K}(x)=H^K-\phi (x), \end{aligned}$$

where $H^K=h(x_{\ell }^{K},y_{\ell }^{K})+\phi (x_{\ell }^{K},y_{\ell }^{K})$ and $\ell =\arg \max _{1\le i\le N_k}\phi (x_{i}^{K},y_{i}^{K})$.

Therefore, we do not require the target still-water equilibria to be known a priori; instead, we can use the source term approximation

$$\begin{aligned} \int _{K}\mathcal {I}\left[ \widehat{\textbf{s}}\right] \cdot \textbf{w}\,dxdy \end{aligned}$$

in the scheme (5.4), where

$$\begin{aligned} \begin{aligned} \widehat{\textbf{s}}=&\,\frac{\textbf{s}(\textbf{U},x,y)}{\textbf{s}(\textbf{u}_e,x,y)}\left( \mathcal {I}[\textbf{f}(\textbf{u}_e)]_x+\mathcal {I}[\textbf{g}(\textbf{u}_e)]_y\right) \\ =&\,\frac{h}{h_{e}^{K}} \begin{pmatrix} 0\\ \mathcal {I}\left[ \frac{1}{2}g(h_{e}^{K})^2\right] _x\\ \mathcal {I}\left[ \frac{1}{2}g(h_{e}^{K})^2\right] _y \end{pmatrix}, \end{aligned} \end{aligned}$$

to preserve the still-water equilibrium exactly.

Remark 5.3

The polytropic equilibrium:

$$\begin{aligned} \frac{p}{\rho ^{\nu }}=\alpha ,\quad \frac{\alpha \nu }{\nu -1}\rho ^{\nu -1}+\phi =\beta ,\quad u=v=0,\quad \nu >1, \alpha , \beta \text { are constants}, \end{aligned}$$

(5.8)

is one of the most basic steady states for the two-dimensional Euler equations:

$$\begin{aligned} {\left\{ \begin{array}{ll} \rho _t+m_x+n_y=0,\\ m_t+\left( \rho u^2+p\right) _x+\left( \rho u v\right) _y=-\rho \phi _x,\\ n_t+\left( \rho u v\right) _x+\left( \rho v^2+p\right) _y=-\rho \phi _y,\\ E_t+\left( (E+p)u\right) _x+\left( (E+p)v\right) _y=-m\phi _x-n\phi _y, \end{array}\right. } \end{aligned}$$

(5.9)

where $\rho$ is the density, $m=\rho u$ and $n=\rho v$ are the momentum in the x and y directions, respectively, (u, v) is the velocity field of the fluid, $E=\frac{1}{2}\rho (u^2+v^2)+\frac{p}{\gamma -1}$ is the total energy, p is the pressure, $\gamma >1$ is the ratio of specific heats, and $\phi$ is the prescribed gravitational potential. For the polytropic equilibrium, the following local reference equilibrium state is available on each cell K:

$$\begin{aligned} \textbf{u}_{e}^{K}= \begin{pmatrix} \rho _{e}^{K}\\ 0\\ 0\\ \frac{p_e^K}{\gamma -1} \end{pmatrix},\quad \rho _{e}^{K} = \left( \frac{\nu -1}{\alpha \nu }(\beta -\phi )\right) ^{\frac{1}{\nu -1}},\quad p_{e}^{K} = \alpha \left( \frac{\nu -1}{\alpha \nu }(\beta -\phi )\right) ^{\frac{\nu }{\nu -1}}, \end{aligned}$$

where $\alpha =\frac{p(x_{\ell }^{K},y_{\ell }^{K})}{\rho (x_{\ell }^{K},y_{\ell }^{K})^\nu }$, $\beta =\frac{\alpha \nu }{\nu -1}\rho (x_{\ell }^{K},y_{\ell }^{K})^{\nu -1}+\phi (x_{\ell }^{K},y_{\ell }^{K})$ and $\ell =\arg \max _{1\le i\le N_k}\phi (x_{i}^{K},y_{i}^{K})$.

Consequently, we do not need the polytropic equilibrium to be known a priori; instead, we can use the source term approximation

$$\begin{aligned} \int _{K}\mathcal {I}\left[ \widehat{\textbf{s}}\right] \cdot \textbf{w}\,dxdy \end{aligned}$$

in the scheme (5.4), where

$$\begin{aligned} \begin{aligned} \widehat{\textbf{s}}=&\frac{\textbf{s}(\textbf{U},x,y)}{\textbf{s}(\textbf{u}_e,x,y)}\left( \mathcal {I}[\textbf{f}(\textbf{u}_e)]_x+\mathcal {I}[\textbf{g}(\textbf{u}_e)]_y\right) \\ =&\begin{pmatrix} 0\\ \frac{\rho }{\rho _{e}^{K}}\mathcal {I}\left[ p_{e}^{K}\right] _x\\ \frac{\rho }{\rho _{e}^{K}}\mathcal {I}\left[ p_{e}^{K}\right] _y\\ -m\phi _x-n\phi _{y} \end{pmatrix}, \end{aligned} \end{aligned}$$

to preserve the polytropic equilibrium exactly.

Note that the above reformulation of the source term was used in [8] on quadrilateral meshes, and its performance has been studied therein. To save space, we will not discuss or test this case further on triangular meshes.

Remark 5.4

Similar to (3.2), the two-dimensional well-balanced scheme (5.4) is designed based on the reformulation of hyperbolic balance laws:

$$\begin{aligned} \textbf{u}_t+\textbf{f}(\textbf{u})_x+\textbf{g}(\textbf{u})_y=\frac{\textbf{s}(\textbf{u},x,y)}{\textbf{s}(\textbf{u}_{e},x,y)}\left( \textbf{f}(\textbf{u}_e)_x+\textbf{g}(\textbf{u}_e)_y\right) , \end{aligned}$$

(5.10)

where $\textbf{u}_e$ is a steady state known a priori.

Alternatively, one can also use the following reformulation [18, 51] to design the well-balanced scheme:

$$\begin{aligned} \textbf{u}_t+\textbf{f}(\textbf{u})_x+\textbf{g}(\textbf{u})_y=\textbf{s}(\textbf{u},x,y)-\textbf{s}(\textbf{u}_{e},x,y)+\left( \textbf{f}(\textbf{u}_e)_x+\textbf{g}(\textbf{u}_e)_y\right) , \end{aligned}$$

resulting in the following implementation:

$$\begin{aligned} \begin{aligned} \int _{K}\textbf{U}_t\cdot \textbf{w}\,dxdy&=\int _{K}\left( \mathcal {I}\left[ \textbf{f}(\textbf{U})\right] \cdot \textbf{w}_x+\mathcal {I}\left[ \textbf{g}(\textbf{U})\right] \cdot \textbf{w}_y\right) \,dxdy\\ &\quad -\int _{\partial K}\mathcal {I}\left[ {\textbf{H}}(\textbf{U}^{\text {int}},\textbf{U}^{\text {ext}},\textbf{n})\right] \cdot \textbf{w}\, ds\\&\quad +\int _{K}\left( \mathcal {I}\left[ \textbf{s}(\textbf{U},x,y)\right] -\mathcal {I}\left[ \textbf{s}(\textbf{u}_e,x,y)\right] +\mathcal {I}[\textbf{f}(\textbf{u}_e)]_x+\mathcal {I}[\textbf{g}(\textbf{u}_e)]_y\right) \\ &\quad \cdot \textbf{w}\,dxdy. \end{aligned} \end{aligned}$$

According to our numerical experiments, this scheme also works well for exactly preserving the target steady states. To save space, we will not showcase the corresponding results.

5.2 Numerical Tests

In this section, we test the performance of the proposed two-dimensional well-balanced DG scheme (5.4) using the classic fourth-order Runge–Kutta time marching method for the shallow water equations (5.2). We take the time step size as $\Delta t=\frac{\text {CFL}}{||\,c+|u|+|v|\, ||_{L^\infty (\Omega )}}h$, where h is the diameter of the maximum inscribed circle of cells, and $\text {CFL}$ constants are 0.2, 0.15, 0.1 for $P^1$, $P^2$ and $P^3$-DG schemes, respectively. Unless otherwise stated, we use the $P^2$-DG scheme and set the gravitational constant $g=9.812$ in the tests. The algorithm is implemented on the Matlab platform in this section. Since it only naturally supports double precision, we only show the error with double precision in the well-balanced tests.

Example 5.1

Accuracy test In this example [52], we test the accuracy of our algorithm. The bottom topography and initial conditions for the shallow water equations are given by:

$$\begin{aligned} \begin{aligned} \phi (x,y)= \sin (2\pi x)+\cos (2\pi y),\quad h(x,y,0) = 10 + \exp \left( \sin (2\pi x)\right) \cos (2\pi y),\\ m(x,y,0) = \sin \left( \cos (2\pi x)\right) \sin (2\pi y),\quad n(x,y,0) = \cos (2\pi x)\cos (\sin (2\pi y)), \end{aligned} \end{aligned}$$

on the computational domain $\Omega =[0,1]^2$ with periodic boundary conditions.

Because a closed-form analytical solution to the problem is not available, we calculate a reference solution using the fifth-order WENO finite difference scheme on a $1280\times 1280$ uniform grid. We compute the well-balanced DG method on uniform triangular grids, which are obtained by splitting uniform $N\times M$ rectangular cells along the diagonal. The reference equilibrium state $\textbf{u}_e$ used in the scheme (5.4) can be chosen arbitrarily. For simplicity, we use the local reference equilibrium of still water, see Remark 5.2.

The $L^1$ errors and orders of convergence at time $T=0.05$ obtained from the DG method are presented in Table 6, from which we can observe high-order accuracy.

Table 6 $L^1$ errors and orders of convergence in the accuracy test 5.1 for the shallow water equations

Example 5.2

Still-water equilibrium In this example [52], we test the well-balanced property of the proposed method for still-water equilibrium. The bottom topography and initial conditions are given by:

$$\begin{aligned} \phi (x,y)=&\,0.8\exp \left( -50\left( (x-0.5)^2+(y-0.5)^2\right) \right) ,\\h(x,y,0)=&1-\phi (x,y),\quad u(x,y,0)=v(x,y,0)=0, \end{aligned}$$

on the computational domain $\Omega =[0,1]^2$ with transmissive boundary conditions.

We compute the well-balanced DG method on uniform triangular grids, which are obtained by splitting uniform $N\times M$ rectangular cells along the diagonal. The surface level $h+\phi$ remains flat in this problem.

The numerical errors at time $T=0.1$ are presented in Table 7, which are clearly at the round-off level.

Table 7 $L^1$ and $L^\infty$ errors on different grids in the well-balanced test 5.2 for a still-water equilibrium at $T=0.1$

Example 5.3

Pseudo-2D moving-water equilibria In this example, we test the well-balanced property of the scheme for pseudo-2D moving-water equilibria. We consider the cases studied in Example 4.2 and Example 4.3, with all variables being constant ($v=0$) along the y direction. The computational domain is $\Omega =[0, 25]\times [0, 10]$ with transmissive boundary conditions. The reference equilibrium is known a priori in these tests. We compute the well-balanced DG method on an unstructured grid, as shown in Fig. 7. The numerical errors at time $T=5$ are presented in Table 8. From the table, we can see the $L^\infty$ errors are at the round-off level. Because the area is large ($|\Omega |=250$), the $L^1$ errors appear somewhat larger.

Fig. 7

Computational grid used in Example 5.3

Table 8 $L^1$ and $L^\infty$ errors for different cases in the well-balanced test 5.3 on an unstructured grid at $T=5$: Case (a) subcritical flow with a continuous bottom; Case (b) supercritical flow with a continuous bottom; Case (c) transcritical flow with a continuous bottom (no shock); Case (d) subcritical flow with a discontinuous bottom; Case (e) supercritical flow with a discontinuous bottom; Case (f) transcritical flow with a discontinuous bottom

Example 5.4

Genuinely 2D moving-water equilibrium In this example, we test the well-balanced property of the proposed scheme for a genuinely two-dimensional moving-water equilibrium. We set the gravitational constant $g=1$ and use a non-flat bottom topography $\phi (x,y)=\frac{1}{2}\cos ^2(x)+\frac{1}{2}\cos ^2(y)$ on the domain $\Omega =[0,2\pi ]^2$ with periodic boundary conditions. One can verify that the following initial conditions

$$\begin{aligned} h(x,y,0)=1,\quad u(x,y,0)=\sin (x)\cos (y),\quad v(x,y,0)=-\cos (x)\sin (y), \end{aligned}$$

constitute a moving-water equilibrium.

We conduct the computation on $20\times 20$, $40\times 40$, and $60\times 60$ uniform triangular grids. The errors at time $T=1.0$ are presented in Table 9, from which we can clearly see that the moving-water equilibrium is preserved exactly up to round-off errors.

Table 9 $L^1$ and $L^\infty$ errors on different grids in the well-balanced test 5.4 at $T=1.0$

Example 5.5

Small perturbation of the genuinely 2D moving-water equilibrium In this example, we test the capability of the well-balanced DG scheme to resolve small perturbations on steady states. We perturb the initial water height with $\delta h(x,y,0)=10^{-3}\exp \left( -10\left( (x-\pi )^2+(y-\pi )^2\right) \right)$ in the moving-water equilibrium presented in Example (5.4), while keeping the other initial conditions and bottom topography the same.

We compute the problem using the proposed well-balanced method, as well as a straightforward discretization of the source term, on three grids with different levels of refinement. The numerical solutions at $T=1$ are presented in Fig. 8. From this comparison, we can clearly observe the excellent resolution of the well-balanced scheme, even on coarse grids.

Fig. 8

Water height computed up to $T=1$ by both well-balanced and non-well-balanced schemes in Example 5.5. Results are obtained from different grids. Row 1: coarse grid; Row 2: medium grid; Row 3: fine grid

6 Concluding Remarks

In this paper, we have established a high order well-balanced discontinuous Galerkin method for the non-hydrostatic equilibria of the shallow water equations and Euler equations. The well-balanced property is achieved through a special reformulation of the source term to mimic the flux gradient. Subsequently, the source term integral is discretized at Gauss-Lobatto points to cancel the flux integral evaluated by the Gauss-Lobatto quadrature and the numerical fluxes. A modified OFDG damping term is added to the scheme to control spurious oscillations near shocks without breaking the well-balanced property. Though the non-hydrostatic equilibria are generally non-polynomials, the steady states are exactly preserved at all Gauss-Lobatto points by the scheme, if the initial conditions are obtained from Lagrange interpolation at Gauss-Lobatto points. As verified by ample numerical tests, the scheme is high-order accurate, well-balanced, and capable of capturing small perturbations of steady states and shocks well.

A The Barycentric Coordinates of the Gauss-Lobatto Type Nodes Adopted on Triangle Cells in Sect. 5

Since the Gauss-Lobatto type nodes coincide with the $(k+1)$ Gauss-Lobatto points on each edge, we list only the barycentric coordinates of the interior points in Table 10.

Table 10 The barycentric coordinates of the Gauss-Lobatto type nodes in the interior of a triangle