# High order well-balanced asymptotic preserving finite difference WENO
schemes for the shallow water equations in all Froude numbers

###### Abstract.

In this paper, high order semi-implicit well-balanced and asymptotic
preserving finite difference WENO schemes are proposed for the shallow water
equations with a non-flat bottom topography. We consider the Froude number
ranging from $\mathcal{O}(1)$ to $0$, which in the zero Froude limit becomes
the “lake equations” for balanced flow without gravity waves. We apply a well-
balanced finite difference WENO reconstruction, coupled with a stiffly
accurate implicit-explicit (IMEX) Runge-Kutta time discretization. The
resulting semi-implicit scheme can be shown to be well-balanced, asymptotic
preserving (AP) and asymptotically accurate (AA) at the same time. Both one-
and two-dimensional numerical results are provided to demonstrate the high
order accuracy, AP property and good performance of the proposed methods in
capturing small perturbations of steady state solutions.

###### Key words and phrases:

shallow water equations; all Froude numbers; finite difference WENO; high
order; asymptotic preserving; well-balanced

Guanlan Huang

School of Mathematical Sciences, Xiamen University

Xiamen, Fujian, 361005, P.R. China

Email<EMAIL_ADDRESS>
Yulong Xing***The work of this author was partially supported by the NSF grant
DMS-1753581.

Department of Mathematics, The Ohio State University

Columbus, OH 43210, USA

Email<EMAIL_ADDRESS>
Tao Xiong†††Corresponding author. The work of this author was partially
supported by NSFC grant No. 11971025, NSF grant of Fujian Province No.
2019J06002, and the Strategic Priority Research Program of Chinese Academy of
Sciences Grant No. XDA25010401.

School of Mathematical Sciences, Xiamen University

Fujian Provincial Key Laboratory of Mathematical Modeling and High-Performance
Scientific Computing

Xiamen, Fujian, 361005, P.R. China

Email<EMAIL_ADDRESS>

##  1\. Introduction

Shallow water equations (SWEs) are widely used in the modeling of water motion
flows in rivers and coastal areas. They have important applications in ocean
currents and hydraulic engineering, see, e.g. [40, 31, 55]. Considering the
water flow in river, reservoir or open channels with a non-flat bottom, the
SWEs can be written as follows:

(1.1) $\left\\{\begin{array}[]{ll}h_{t}+\nabla\cdot(h{\bf
u})=0,\\\\[8.53581pt] (h{\bf u})_{t}+\nabla\cdot(h{\bf u}\otimes{\bf
u})+g\nabla(h^{2}/2)=-gh\nabla b,\end{array}\right.$

where $h$ is the depth of the water layer, ${\bf u}$ is the flow velocity,
defined on a time-space domain $(t,{\bf x})\in\mathbb{R}^{+}\times\Omega$. $g$
is the gravitational constant and $b({\bf x})$ is the bottom topography which
is independent of time. $\otimes$ denotes the Kronecker product. When the
bottom is flat, this system is equivalent to the isentropic Euler system in
the homogeneous case. However, the geometrical source term changes the
property of the system when a non-flat bottom topography is taken into
account.

Many shock capturing schemes with explicit time discretizations have been
developed to solve the SWEs with source term (1.1), including high order
finite difference [58, 62, 20, 36, 59], finite volume [1, 43, 10, 34, 44, 64,
38], residual distribution methods [47, 46] and discontinuous Galerkin schemes
[65, 61, 60, 67], and many references therein. When solving the SWEs with
source term numerically, it is important to preserve the exact conservation
property (C-property) [35], namely, the nonzero flux gradient should be
exactly balanced by the source term in the case of a stationary water. Such
schemes are named well-balanced methods. During the past few decades, there
have been extensive studies on the design and analysis of well-balanced
methods for various hyperbolic equations with source terms. For the SWEs
(1.1), the still-water stationary solution takes the form

(1.2) $h{\bf u}=0,\quad h+b=\text{Const.}$

Well-balanced schemes for the SWEs are able to capture small perturbations of
the hydrostatic or nearly hydrostatic flows on a coarse mesh, and we refer to
the above list of literatures and the survey papers [63, 33] for more
discussions.

On the other hand, if we choose a characteristic length $l_{0}$, a
characteristic depth $h_{0}$, a characteristic velocity $U_{0}$ and a
characteristic time $t_{0}$, we can define the following dimensionless
variables

(1.3) $\hat{{{\bf x}}}=\frac{{{\bf
x}}}{\ell_{0}},\quad\hat{h}=\frac{h}{h_{0}},\quad\hat{{\bf u}}=\frac{{\bf
u}}{U_{0}},\quad\hat{t}=\frac{t}{t_{0}},\quad\hat{b}=\frac{b}{h_{0}},$

with which, the SWEs (1.1) can be non-dimensionalized as follows:

(1.4) $\left\\{\begin{array}[]{ll}{\text{Sr}}\,h_{t}+\nabla\cdot(h{\bf
u})=0,\\\\[8.53581pt] {\text{Sr}}\,(h{\bf u})_{t}+\nabla\cdot(h{\bf
u}\otimes{\bf
u})+\frac{1}{{\text{Fr}}^{2}}\nabla(h^{2}/2)=-\frac{1}{{\text{Fr}}^{2}}h\nabla
b,\end{array}\right.$

where we drop the hat of the dimensionless variables for ease of presentation.
The Strouhal number Sr and the Froude number Fr are defined as

(1.5)
${\text{Sr}}:=\frac{\ell_{0}}{t_{0}U_{0}},\quad{\text{Fr}}:=\frac{U_{0}}{\sqrt{gh_{0}}}.$

In case of low Froude number flows, for which the flow velocities are
systematically small as compared to the velocity of gravity waves, a reference
asymptotic expansion parameter $\varepsilon$ can be introduced according to
the Froude number, via

(1.6) ${\text{Fr}}=\varepsilon^{\alpha}\ll 1,$

with $\alpha$ chosen depending on the particular flow regime to be considered
[32].

In this work, we focus on flows over advective time scale where
${\text{Sr}}=1$, and assume ${\text{Fr}}=\varepsilon$, namely $\alpha=1$ for
the inviscid balanced flow over the topography, so that the dimensionless
equations (1.4) become

(1.7) $\left\\{\begin{array}[]{ll}h_{t}+\nabla\cdot(h{\bf
u})=0,\\\\[8.53581pt] (h{\bf u})_{t}+\nabla\cdot(h{\bf u}\otimes{\bf
u})+\frac{1}{\varepsilon^{2}}\nabla(h^{2}/2)=-\frac{1}{\varepsilon^{2}}h\nabla
b.\end{array}\right.$

The system is hyperbolic, and its eigenvalues in the direction ${\bf{n}}$ are
$\lambda_{1}={\bf u}\cdot{\bf{n}}+c/\varepsilon$ and $\lambda_{2}={\bf
u}\cdot{\bf{n}}-c/\varepsilon$, with $c=\sqrt{h}$ being the scaled speed of
sound.

One could directly apply the well-balanced shock capturing schemes to the
dimensionless system (1.7), however, due to the fact that the characteristic
speed $\lambda_{1,2}$ is inversely proportional to the Froude number
$\varepsilon$, the time step constraint of an explicit time discretization
satisfies

$\Delta t=\text{CFL}\frac{\Delta x}{\max(|{\bf
u}|+c/\varepsilon)}\sim\varepsilon\Delta x,$

where $\Delta t$ is the time step size, $\Delta x$ is the mesh size and CFL is
the time stability CFL number. As the Froude number $\varepsilon$ approaching
to zero, this leads to the stiffness in time, which is the same as the low
Mach flows, see e.g. [24, 14, 12]. For low Mach flows, preconditioning
techniques are usually applied to release the small time step condition and
cure large numerical viscosities in the shock capturing schemes [54, 57, 15,
41, 11]. Such techniques, however, are effectively applicable only if the Mach
numbers are not too small. On the other hand, naive implicit time
discretizations of these shock capturing schemes result in the fully nonlinear
systems, which are very inefficient to solve and sometimes may not be able to
converge to the correct asymptotic limit.

In between, many semi-implicit schemes are developed, e.g., for low Mach (all
Mach) Euler and Navier-Stokes equations [14, 24, 50, 18, 7, 5, 12, 51, 17, 66,
16, 52, 8, 6], and for low Froude shallow water equations [22, 21, 53, 2, 56,
37], and many references therein. Among them, one type of methods which can
ensure the correct asymptotic limit is the asymptotic preserving (AP) scheme.
AP schemes were originally introduced in [27] for multiscale kinetic
equations, namely, the discretized scheme for a stiff PDE can converge to a
consistent discretization of its limiting equation, under unresolved time step
and mesh sizes, with uniform stability. For a recent review of AP schemes and
their applications, see [25]. For Euler or Navier-Stokes equations will all-
Mach number, AP schemes can well capture its corresponding incompressible
limit as the Mach number approaching to zero [14, 24, 12, 42, 3, 5]. AP
schemes have also been applied to the shallow water equations with the low
Froude number limit [2, 19, 13, 39, 37]. However, for the shallow water
equations with an irregular bottom topography, most of current schemes either
designed focus only in the low Froude number regime, or with up to second
order accuracy.

In this paper, we propose high order well-balanced asymptotic preserving
weighted essentially non-oscillatory (WENO) schemes for the shallow water
equations with a non-flat bottom topography and all Froude numbers. For an
irregular bottom with $b({\bf x})\neq 0$, the nonzero flux gradient and the
source term are both scaled by the Froude number in (1.7). Here, considering
the still water equilibrium (1.2), it is important to preserve both the well-
balanced property and the low Froude limit for the dimensionless system (1.7).
A close work along this line is the second order well-balanced asymptotic
preserving scheme developed and carefully analyzed by Liu in [37], which is
based on a split system for the pre-balanced shallow water equations,
following the idea for all Mach flows in [24]. In our work, we will combine
the high order AP schemes developed for the isentropic Euler and full Euler
systems with all Mach numbers [5, 6] in the spirit of [14, 50], with the well
balanced finite difference WENO schemes [62], to achieve a high order well-
balanced asymptotic preserving scheme for the shallow water equations with a
source term. We start by constructing a first order semi-implicit scheme.
Similar to the hydrostatic pressure $p_{2}$ introduced for the pressure in the
all-Mach flow [5, 6], here an $H_{2}$ term corresponding to the variation from
a constant water surface level with respect to the total water surface $H=h+b$
is introduced. We first solve $H_{2}$ from an elliptic (or Helmholtz)
equation, which is formed from a semi-implicit time discretization. After
evaluating $H_{2}$, we can update the momentum $h{\bf u}$ and then $h$. In
this way, by utilizing a well balanced flux reconstruction in the updating of
$h$ and $h{\bf u}$, we can show that our first order semi-discrete scheme
achieves the well balanced and AP properties at the same time. With the aid of
a multi-stage explicit-implicit (IMEX) Runge-Kutta time discretization for a
partitioned autonomous system, high order semi-implicit schemes can be
obtained. Corresponding suitable high order spatial discretizations can also
be constructed. Specifically high order well-balanced finite difference WENO
reconstruction [62] for convection terms are used in this paper, with high
order central difference discretizations of second order and mixed derivatives
in the elliptic (or Helmholtz) equation of $H_{2}$. The resulting high order
semi-implicit scheme is showed to satisfy the well-balanced, AP and
asymptotically accurate (AA) properties simultaneously, namely, the scheme is
not only consistent (AP property) but also preserves the order of accuracy in
time (AA property) in the stiff limit as $\varepsilon\rightarrow 0$ [45].

The rest of the paper is organized as follows. In Section 2, the low Froude
limit of the SWEs is revisited. In Section 3, a well balanced AP scheme based
on a first order semi-implicit scheme is first described and then generalized
to high order methods. The analysis of well-balanced property, as well as AP
and AA properties, follows afterward. Numerical experiments are presented in
Section 4, which demonstrate the good performance of the high order well-
balanced AP scheme in nearly hydrostatic flows and for a range of the Froude
numbers including the zero Froude number limit. Conclusions are made in
Section 5.

##  2\. Low Froude number limit for SWEs

Let us denote $H=h+b$ as the water surface level, and the system (1.7) can be
written as

(2.1) $\left\\{\begin{array}[]{ll}h_{t}+\nabla\cdot(h{\bf
u})=0,\\\\[8.53581pt] (h{\bf u})_{t}+\nabla\cdot(h{\bf u}\otimes{\bf
u})+\frac{1}{\varepsilon^{2}}h\nabla H=0.\end{array}\right.$

We start with the following single-scale expansions of the solutions $h$ and
${\bf u}$, in terms of $\varepsilon$,

(2.2) $\left\\{\begin{array}[]{ll}h({\bf x},t)=h_{0}({\bf x},t)+\varepsilon
h_{1}({\bf x},t)+\varepsilon^{2}h_{2}({\bf x},t)+\cdots,\\\\[8.53581pt] {\bf
u}({\bf x},t)={\bf u}_{0}({\bf x},t)+\varepsilon{\bf u}_{1}({\bf
x},t)+\varepsilon^{2}{\bf u}_{2}({\bf x},t)+\cdots.\end{array}\right.$

Since $H=h+b$ with $b=b({\bf x})$ being time independent, we have

(2.3) $H({\bf x},t)=h_{0}({\bf x},t)+b({\bf x})+\varepsilon h_{1}({\bf
x},t)+\varepsilon^{2}h_{2}({\bf x},t)+\cdots$

Substituting (2.2) and (2.3) into (2.1), equating to zero for different orders
of $\varepsilon$, we have

  * •

$\mathcal{O}(\varepsilon^{-2})$

(2.4) $h_{0}\nabla(h_{0}+b)=0,$

  * •

$\mathcal{O}(\varepsilon^{-1})$

(2.5) $h_{1}\nabla(h_{0}+b)+h_{0}\nabla h_{1}=0,$

  * •

$\mathcal{O}(\varepsilon^{0})$

(2.6) $\left\\{\begin{array}[]{ll}(h_{0})_{t}+\nabla\cdot(h_{0}{\bf
u}_{0})=0,\\\\[8.53581pt] (h_{0}{\bf u}_{0})_{t}+\nabla\cdot(h_{0}{\bf
u}_{0}\otimes{\bf u}_{0})+h_{2}\nabla(h_{0}+b)+h_{1}\nabla h_{1}+h_{0}\nabla
h_{2}=0.\end{array}\right.$

Here for simplicity, no dry area is considered to exist in the domain so that
$h_{0}\neq 0$. Therefore, from (2.4) we obtain

(2.7) $h_{0}+b\equiv H_{0}(t),$

namely, $h_{0}+b$ is constant in space. It follows from (2.5) that
$h_{1}\equiv H_{1}(t)$ is also constant in space. Since the bottom topography
$b$ is assumed to be time independent, from (2.6), we have

(2.8a) $\nabla\cdot(h_{0}{\bf u}_{0})=-\frac{dH_{0}(t)}{dt},$ (2.8b)
$(h_{0}{\bf u}_{0})_{t}+\nabla\cdot(h_{0}{\bf u}_{0}\otimes{\bf
u}_{0})+h_{0}\nabla h_{2}=0.$

Now integrating the equation (2.8a) over the spatial domain $\Omega$, it
yields

(2.9)
$\frac{dH_{0}(t)}{dt}=-\frac{1}{|\Omega|}\int_{\Omega}\nabla\cdot(h_{0}{\bf
u}_{0})d\sigma=-\frac{1}{|\Omega|}\int_{\partial\Omega}h_{0}{\bf
u}_{0}\cdot{\bf n}ds,$

where ${\bf n}$ is the unit outward normal vector along $\partial\Omega$,
namely the time change of the total water height is given by the total flux of
water across the domain boundary. (2.8) and (2.9) form the classical zero
Froude number shallow water equations, also known as the “lake equations” [23,
32]. If considering the no-slip ${\bf u}\cdot{\bf n}=0$ or periodic boundary
conditions, we further get $\int_{\Omega}\nabla\cdot(h_{0}{\bf u}_{0})\,d{\bf
x}=\int_{\partial\Omega}h_{0}{\bf u}_{0}\cdot{\bf n}\,ds=0$. This implies
$H_{0}$ is constant both in space and time, i.e. $H_{0}=$ Const. The same
conclusion can also be derived for $H_{1}$. Therefore, the “lake equations”
further reduce to:

(2.10) $\left\\{\begin{array}[]{l}\nabla\cdot(h_{0}{\bf u}_{0})=0,\quad
h_{0}+b=H_{0}=\text{Const.}\\\\[8.53581pt] \displaystyle\partial_{t}(h_{0}{\bf
u}_{0})+\nabla\cdot\left(h_{0}{\bf u}_{0}\otimes{\bf u}_{0}\right)+h_{0}\nabla
h_{2}=0.\end{array}\right.$

A rigrous convergence analysis for the zero Froude limit from (2.1) to (2.10)
is very demanding, and we refer to [28, 29] for such a rigorous study in the
low Mach limit.

##  3\. Numerical schemes

In this section, we will construct and analyze a class of high order finite
difference schemes with the AP and well-balanced properties for the shallow
water equations (2.1) with a range of Froude numbers. The SWEs in the form of
(2.1) are very close to the isentropic Euler equations with all-Mach number,
see e.g. [14, 5]. However, it differs in the zero Froude limit, where in the
isentropic Euler system, $\rho_{0}$ (corresponding to $h_{0}$ here) is
constant, and it has the divergence free velocity field $\nabla\cdot{\bf
u}_{0}=0$. Here, $h_{0}$ is not a constant directly due to the appearance of
source term. For the isentropic Euler equation in the zero Mach limit, the
hydrodynamic pressure $p_{2}$ (corresponding to $H_{2}$ here) plays a role as
a Lagrangian multiplier to ensure the divergence free condition, fortunately
$H_{2}$ performs similarly in this setting. In the following, we will extend
the high order semi-implicit finite difference WENO schemes developed in [5]
to solve (2.1). We will analyze that the scheme can capture the zero Froude
number shallow water equations, or the lake equations (2.10) with no-slip or
periodic boundary conditions, namely the scheme is asymptotic preserving.

For the shallow water equations with non-flat bottom topography, the well-
balanced property is another important one, especially for capturing small
perturbations of a still water equilibrium [62]. We will adopt the well-
balanced finite difference WENO reconstruction technique as developed in [62],
tailored to our semi-implicit time discretization. We will show that under our
semi-implicit framework, the well-balanced property can also be obtained.

###  3.1. First order semi-implicit scheme

We start with presenting a first order semi-implicit time discretization,
while keeping space continuous at this moment. The first order semi-implicit
implicit-explicit (IMEX) scheme for (2.1) is given as follows

(3.1) $\left\\{\begin{aligned} &\frac{h^{n+1}-h^{n}}{\Delta
t}+\nabla\cdot(h{\bf u})^{n+1}=0,\\\ &\frac{(h{\bf u})^{n+1}-(h{\bf
u})^{n}}{\Delta t}+\nabla\cdot\left(\frac{h{\bf u}\otimes h{\bf
u}}{h}\right)^{n}+\frac{1}{\varepsilon^{2}}h^{n+1}\nabla
H^{n+1}=0.\end{aligned}\right.$

Notice that $H^{n+1}=h^{n+1}+b$. From the second equation of (3.1), we can
first express $(h{\bf u})^{n+1}$ in terms of $h^{n+1}$, $H^{n+1}$ and other
variables at time level $t^{n}$. Substituting it into the first equation of
(3.1), we get

(3.2) $\left\\{\begin{aligned} &\frac{h^{n+1}-h^{n}}{\Delta
t}+\nabla\cdot(h{\bf u})^{n}-\Delta t\nabla^{2}:\left(\frac{h{\bf u}\otimes
h{\bf u}}{h}\right)^{n}-\frac{\Delta
t}{\varepsilon^{2}}\nabla\cdot(h^{n+1}\nabla H^{n+1})=0,\\\ &\frac{(h{\bf
u})^{n+1}-(h{\bf u})^{n}}{\Delta t}+\nabla\cdot\left(\frac{h{\bf u}\otimes
h{\bf u}}{h}\right)^{n}+\frac{1}{\varepsilon^{2}}h^{n+1}\nabla
H^{n+1}=0,\end{aligned}\right.$

where $:$ is the tensor double dot product. The first equation of (3.2)
appears to be a nonlinear system for $h^{n+1}$, as $H^{n+1}=h^{n+1}+b$. To
avoid it, a slight modification of the first equation yields

(3.3) $\frac{h^{n+1}-h^{n}}{\Delta t}+\nabla\cdot(h{\bf u})^{n}-\Delta
t\nabla^{2}:\left(\frac{h{\bf u}\otimes h{\bf u}}{h}\right)^{n}-\frac{\Delta
t}{\varepsilon^{2}}\nabla\cdot(h^{n}\nabla H^{n+1})=0,$

which is now a linear equation for the unknown function $h^{n+1}$. Similarly,
the corresponding term in the second equation of (3.2) can be replaced by
$\frac{1}{\varepsilon^{2}}h^{n}\nabla H^{n+1}$, which is easier for the
generalization to high order described in the following subsection.

To deal with the stiff diffusive term
$\frac{1}{\varepsilon^{2}}\nabla\cdot(h^{n}\nabla H^{n+1})$, in this work, we
consider the no-slip or periodic boundary conditions, namely, $h_{0}+b=H_{0}$
and $h_{1}=H_{1}$ are both constants in the asymptotic expansion (2.3). We may
now introduce a water surface perturbation $H_{2}$, defined as

(3.4) $H_{2}=\frac{H-\bar{H}}{\varepsilon^{2}},$

where $\bar{H}$ denotes the spatial average of the water surface level $H$
(computed from $h+b$). In this way, the term $H_{2}$ in the zero Froude
shallow water limit converges to $h_{2}$, which remains finite. Numerically,
we take $\bar{H}$ as the spatial average of $H^{n}$, that is

(3.5) $H^{n+1}=\bar{H}^{n}+\varepsilon^{2}H^{n+1}_{2}\quad\text{ and }\quad
h^{n+1}=\bar{H}^{n}-b+\varepsilon^{2}H^{n+1}_{2},$

so that we obtain a linear elliptic equation for $H^{n+1}_{2}$ from (3.3)

(3.6) $\varepsilon^{2}H^{n+1}_{2}-\Delta t^{2}\nabla\cdot(h^{n}\nabla
H^{n+1}_{2})=h^{*},$

with

(3.7) $h^{*}=H^{n}-\bar{H}^{n}-\Delta t\left(\nabla\cdot(h{\bf u})^{n}-\Delta
t\nabla^{2}:\left(\frac{h{\bf u}\otimes h{\bf u}}{h}\right)^{n}\right).$

After obtaining $H^{n+1}_{2}$ from (3.6), $h^{n+1}$ and $H^{n+1}$ can be
updated from (3.5). In equation (3.2), we can replace
$1/\varepsilon^{2}h^{n+1}\nabla H^{n+1}$ by $h^{n}\nabla H^{n+1}_{2}$, which
leads to

(3.8) $\left\\{\begin{aligned} &\frac{h^{n+1}-h^{n}}{\Delta
t}+\nabla\cdot(h{\bf u})^{n}-\Delta t\nabla^{2}:\left(\frac{h{\bf u}\otimes
h{\bf u}}{h}\right)^{n}-\nabla\cdot(h^{n}\nabla H_{2}^{n+1})=0,\\\
&\frac{(h{\bf u})^{n+1}-(h{\bf u})^{n}}{\Delta t}+\nabla\cdot\left(\frac{h{\bf
u}\otimes h{\bf u}}{h}\right)^{n}+h^{n}\nabla
H_{2}^{n+1}=0.\end{aligned}\right.$

We can solve the second equation of (3.8) for $(h{\bf u})^{n+1}$. This semi-
implicit treatment in time can ensure the right asymptotic limit as the Froude
number $\varepsilon\rightarrow 0$, which is known as the AP property. We will
analyze it afterward. In general, direct computing $h^{n+1}$ from (3.5) cannot
preserve exact mass conservation, and we may further update $h^{n+1}$ using
the first equation of (3.1) with the available $(h{\bf u})^{n+1}$.

Next we will discuss the spatial discretizations according to the first order
semi-implicit time discretization. The main guidance is to preserve the
equilibrium state for a still water when $H=h+b=$ Const. and $h{\bf u}={\bf
0}$, and also avoid excessive numerical viscosity inversely proportional to
the Froude number $\varepsilon$. We follow both the well-balanced finite
difference scheme developed in [62], and the spatial discretizations for the
all-Mach isentropic Euler equations in [5]. First or second low order
discretizations will be described first, and high order extensions will be
presented afterward.

To preserve the still water equilibrium for the water surface level $H$, it is
more convenient to rewrite the first equation of (3.1) in a pre-balanced form,
namely

(3.9) $\frac{H^{n+1}-H^{n}}{\Delta t}+\nabla\cdot(h{\bf u})^{n+1}=0,$

which is equivalent to the original equation since $b$ is independent of time.
In the case of still-water equilibrium (1.2), to preserve the water surface
level $H=$ Const, it requires that no numerical viscosity should be presented
in the numerical approximation of the flux term $\nabla\cdot(h{\bf u})^{n+1}$.
Therefore, the numerical viscosity term should depend on $H$ instead of $h$,
and a Lax-Friedrichs flux for $\nabla\cdot(h{\bf u})$ (we drop the superindex
$n+1$ for brevity) is defined as follows:

(3.10a)
$(\widehat{hu})_{i+\frac{1}{2},j}=\frac{1}{2}\left[(hu)_{i+1,j}+(hu)_{i,j}-\alpha^{x}_{i,j}(H_{i+1,j}-H_{i,j})\right],$
(3.10b)
$(\widehat{hv})_{i,j+\frac{1}{2}}=\frac{1}{2}\left[(hv)_{i,j+1}+(hv)_{i,j}-\alpha^{y}_{i,j}(H_{i,j+1}-H_{i,j})\right],$

where $\alpha^{x}_{i,j}=\max_{h,hu}(|u|+\min(1,1/\varepsilon)\sqrt{h})$ and
$\alpha^{y}_{i,j}=\max_{h,hv}(|v|+\min(1,1/\varepsilon)\sqrt{h})$ are the
local viscosity coefficients along $x$ and $y$ respectively. We denote

(3.11) $\nabla_{LF}\cdot(h{\bf
u})=\frac{(\widehat{hu})_{i+\frac{1}{2},j}-(\widehat{hu})_{i-\frac{1}{2},j}}{\Delta
x}+\frac{(\widehat{hv})_{i,j+\frac{1}{2}}-(\widehat{hv})_{i,j-\frac{1}{2}}}{\Delta
y}.$

It is easy to see $\nabla_{LF}\cdot(h{\bf u})=0$ if $H=h+b=$ Const. and ${\bf
u}=(u,v)=(0,0)$.

For the second equation of (3.8), it does not matter very much how to choose
the numerical fluxes since the numerical viscosity depends on $h{\bf u}$ which
is $0$ for still water, e.g., a local Lax-Friedrichs flux for the second term
$\nabla\cdot\left(h{\bf u}\otimes h{\bf u}/h\right)$ is

(3.12a)
$\left(\widehat{hu^{2}}\right)_{i+\frac{1}{2},j}=\frac{1}{2}\left[\left(hu^{2}\right)_{i+1,j}+\left(hu^{2}\right)_{i,j}-\alpha^{x}_{i,j}\Big{(}(hu)_{i+1,j}-(hu)_{i,j}\Big{)}\right],$
(3.12b)
$\left(\widehat{huv}\right)_{i,j+\frac{1}{2}}=\frac{1}{2}\left[\left(huv\right)_{i,j+1}+\left(huv\right)_{i,j}-\alpha^{y}_{i,j}\Big{(}(hu)_{i,j+1}-(hu)_{i,j}\Big{)}\right],$
(3.12c)
$\left(\widehat{huv}\right)_{i+\frac{1}{2},j}=\frac{1}{2}\left[\left(huv\right)_{i+1,j}+\left(huv\right)_{i,j}-\alpha^{x}_{i,j}\Big{(}(hv)_{i+1,j}-(hv)_{i,j}\Big{)}\right],$
(3.12d)
$\left(\widehat{hv^{2}}\right)_{i,j+\frac{1}{2}}=\frac{1}{2}\left[\left(hv^{2}\right)_{i,j+1}+\left(hv^{2}\right)_{i,j}-\alpha^{y}_{i,j}\Big{(}(hv)_{i,j+1}-(hv)_{i,j}\Big{)}\right],$

and $\alpha^{x}_{i,j}$ and $\alpha^{y}_{i,j}$ are the local viscosity
coefficients which can be taken the same as above.

Notice that preserving the still water equilibrium in (3.8) is to require
$(h{\bf u})^{n+1}={\bf 0}$, which can be satisfied from requiring $h^{n}\nabla
H^{n+1}_{2}=0$. However, a straightforward numerical discretization of the
term $h\nabla H_{2}$ may lead to a nonconservative discretization, even in the
special case of $b=0$ when such term should be treated in the conservative
manner. This issue has been addressed in the well-balanced WENO methods
studied in [62]. By adopting such idea to decompose the source term, at the
continuous level using the relation (3.4) we can rewrite this term as

(3.13) $h\nabla H_{2}=\frac{1}{\varepsilon^{2}}h\nabla
H=\frac{1}{\varepsilon^{2}}\nabla\left(\frac{1}{2}h^{2}-\frac{1}{2}b^{2}\right)+\frac{1}{\varepsilon^{2}}H\nabla
b=\nabla\left(\bar{H}H_{2}+\frac{1}{2}\varepsilon^{2}H^{2}_{2}-H_{2}b\right)+H_{2}\nabla
b.$

Taking $H_{2}$ as $H^{n+1}_{2}$ and using central differences for a low order
spatial discretization for both terms, that is

$\nabla_{C}\bigg{(}\bar{H}^{n}H^{n+1}_{2}+\frac{1}{2}\varepsilon^{2}(H^{n+1}_{2})^{2}-H^{n+1}_{2}b\bigg{)}+H^{n+1}_{2}\nabla_{C}b,$

which is still $0$ when $H_{2}^{n+1}\equiv 0$.

The remaining spatial discretizations for (3.6) and (3.7) are as follows. We
use central difference discretization for the second order derivatives terms,
denoted with subindex $C$, where

(3.14) $\displaystyle\nabla^{2}_{C}:\left(\frac{h{\bf u}\otimes h{\bf
u}}{h}\right)=$ $\displaystyle\frac{1}{\Delta
x^{2}}\Big{[}(hu^{2})_{i+1,j}-2(hu^{2})_{i,j}+(hu^{2})_{i-1,j}\Big{]}$
$\displaystyle+\frac{1}{2\Delta x\Delta
y}\Big{[}\Big{(}(huv)_{i+1,j+1}-(huv)_{i-1.j+1}\Big{)}-\Big{(}(huv)_{i+1,j-1}-(huv)_{i-1,j-1}\Big{)}\Big{]}$
$\displaystyle+\frac{1}{\Delta
y^{2}}\Big{[}(hv^{2})_{i,j+1}-2(hv^{2})_{i,j}+(hv^{2})_{i,j-1}\Big{]}.$

For the term $\nabla\cdot(h^{n}\nabla H^{n+1})$ expressing in the form

$\nabla\cdot(h^{n}\nabla
H^{n+1})=\partial_{x}(h^{n}\partial_{x}H^{n+1})+\partial_{y}(h^{n}\partial_{y}H^{n+1}),$

we may take a compact central difference for terms like $(a(x,y)q_{x})_{x}$ at
the grid point $(x_{i},y_{j})$

$(a(x,y)q_{x})_{x}\Big{|}_{(x_{i},y_{j})}=\frac{1}{\Delta
x^{2}}(a_{i-1,j},a_{i,j},a_{i+1,j})\begin{pmatrix}\frac{1}{2}&-\frac{1}{2}&0\\\
\frac{1}{2}&-1&\frac{1}{2}\\\
0&-\frac{1}{2}&\frac{1}{2}\end{pmatrix}\begin{pmatrix}q_{i-1,j}\\\ q_{i,j}\\\
q_{i+1,j}\end{pmatrix},$

and similar approximation can be done for $(a(x,y)q_{y})_{y}$ along the $y$
direction. This will form a positive definite linear system for the left side
of (3.6), if $h^{n}$ keeps positive. We denote the numerical approximation of
$\nabla\cdot(h^{n}\nabla H^{n+1})$ by $\nabla_{C^{2}}\cdot(h^{n}\nabla
H^{n+1})$. Lastly $\nabla\cdot(h{\bf u})^{n}$ is discretized the same as in
(3.11), so that on the right side of (3.6)

(3.15) $h^{*}=H^{n}-\bar{H}^{n}-\Delta t\left(\nabla_{LF}\cdot(h{\bf
u})^{n}-\Delta t\nabla_{C}^{2}:\left(\frac{h{\bf u}\otimes h{\bf
u}}{h}\right)^{n}\right),$

which is clearly $0$ for $H=$ Const. and $h{\bf u}={\bf 0}$. With such
discretizations, solving $H^{n+1}_{2}$ from (3.6) yields $H^{n+1}_{2}\equiv
0$, so the well-balanced property for the still water is well preserved.

We now summarize the first order semi-implicit scheme as follows:

(3.16) $\left\\{\begin{aligned} &\varepsilon^{2}H^{n+1}_{2}-\Delta
t^{2}\nabla_{C^{2}}\cdot(h^{n}\nabla H^{n+1}_{2})=H^{n}-\bar{H}^{n}-\Delta
t\left(\nabla_{LF}\cdot(h{\bf u})^{n}-\Delta t\nabla^{2}_{C}:\left(\frac{h{\bf
u}\otimes h{\bf u}}{h}\right)^{n}\right),\\\ &\frac{(h{\bf u})^{n+1}-(h{\bf
u})^{n}}{\Delta t}+\nabla_{LF}\cdot\left(\frac{h{\bf u}\otimes h{\bf
u}}{h}\right)^{n}+\nabla_{C}\bigg{(}\bar{H}^{n}H^{n+1}_{2}+\frac{1}{2}\varepsilon^{2}(H^{n+1}_{2})^{2}-H^{n+1}_{2}b\bigg{)}+H^{n+1}_{2}\nabla_{C}b=0,\\\
&\frac{h^{n+1}-h^{n}}{\Delta t}+\nabla_{LF}\cdot(h{\bf
u})^{n+1}=0,\end{aligned}\right.$

which is performed in a sequential way.

###  3.2. High order semi-implicit scheme

To extend the first order semi-implicit scheme to high order, we follow a
similar procedure as described in [5, 7]. For ease of presentation, we keep
space continuous first. Let’s write (2.1) as an autonomous system

(3.17) $U_{t}=\mathcal{H}(U,U),$

where $U=(h,h{\bf u})^{T}$ and
$\mathcal{H}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a
sufficiently regular mapping. We use two different arguments for $U$ with
different treatments, one is explicit with subindex “E” and the other is
implicit with subindex “I”, that is $U_{E}=(h_{E},(h{\bf u})_{E})^{T}$ and
$U_{I}=(h_{I},(h{\bf u})_{I})^{T}$, and we solve

(3.18)
$\left\\{\begin{array}[]{l}U^{\prime}_{E}=\mathcal{H}(U_{E},U_{I}),\\\\[8.53581pt]
U^{\prime}_{I}=\mathcal{H}(U_{E},U_{I}),\end{array}\right.$

where we define

(3.19) $\mathcal{H}(U_{E},U_{I})=\begin{pmatrix}-\nabla\cdot(h{\bf u})_{I},\\\
\,\\\ -\nabla\cdot\left(h{\bf u}\otimes{\bf u}\right)_{E}-h_{E}\nabla
H_{I,2}\end{pmatrix}.$

$H_{I,2}$ is defined similarly as in (3.4)

(3.20)
$H_{I,2}=\frac{H_{I}-\bar{H}_{E}}{\varepsilon^{2}}=\frac{h_{I}+b-\bar{H}_{E}}{\varepsilon^{2}},$

and $\bar{H}_{E}$ is the spatial average of $h_{E}+b$. For the first order
semi-implicit scheme, $U_{E}=U^{n}=(h^{n},(h{\bf u})^{n})^{T}$ and
$U_{I}=U^{n+1}=(h^{n+1},(h{\bf u})^{n+1})^{T}$.

For the partitioned system (3.18), we need to apply an IMEX Runge-Kutta time
discretization with a double Butcher $tableau$ [9],

(3.21) $\begin{array}[]{c|c}\tilde{c}&\tilde{A}\\\
\hline\cr\vspace{-0.25cm}\hfil\\\ &\tilde{b^{T}}\end{array}\ \ \ \ \
\qquad\begin{array}[]{c|c}{c}&{A}\\\ \hline\cr\vspace{-0.25cm}\hfil\\\
&{b^{T}}\end{array},$

where $\tilde{A}=(\tilde{a}_{ij})$ is an $s\times s$ matrix for an explicit
scheme, with $\tilde{a}_{ij}=0$ for $j\geq i$ and $A=({a}_{ij})$ is an
$s\times s$ matrix for an implicit scheme. For the implicit part of the
methods, we use a diagonally implicit scheme, i.e. $a_{ij}=0$, for $j>i$, in
order to guarantee simplicity and efficiency in solving the algebraic
equations corresponding to the implicit part of the discretization. The
vectors $\tilde{c}=(\tilde{c}_{1},...,\tilde{c}_{s})^{T}$,
$\tilde{b}=(\tilde{b}_{1},...,\tilde{b}_{s})^{T}$, and
$c=(c_{1},...,c_{s})^{T}$, $b=(b_{1},...,b_{s})^{T}$ complete the
characterization of the scheme. The coefficients $\tilde{c}$ and $c$ are given
by the usual relation

(3.22) $\displaystyle\tilde{c}_{i}=\sum_{j=1}^{i-1}\tilde{a}_{ij},\ \ \
c_{i}=\sum_{j=1}^{i}a_{ij}.$

For the first order semi-implicit scheme, it corresponds to $s=1$, and the
double Butcher Tableau is

$\begin{array}[]{c|c}0&0\\\
\hline\cr&1\end{array}\qquad\qquad\begin{array}[]{c|c}1&1\\\
\hline\cr&1\end{array},$

namely $U_{E}=U^{n}$ and $U_{I}=U^{n+1}$.

For a high order semi-implicit scheme, a multi-stage IMEX Runge-Kutta is
needed, usually it is characterized as the triplet $(s,\sigma,p)$, for the
number of stages of the implicit scheme ($s$), the number of stages of the
explicit scheme ($\sigma$) and the order of the scheme ($p$). Here we adopt
the IMEX scheme as constructed in [6], which we require $\sigma=s$ with $s$
stages for both implicit and explicit parts, and $\tilde{c}_{i}=c_{i}$ for
$i=2,\cdots,s$.

Now, we may update the solutions as follows. Starting from
$U_{E}^{(0)}=U_{I}^{(0)}=U^{n}$, for inner stages $i=1\text{ to }s$:

  * •

First update the solution $U_{E}^{(i)}$ for the explicit part

(3.23) $U_{E}^{(i)}=U^{n}+\Delta
t\sum^{i-1}_{j=1}\tilde{a}_{ij}\mathcal{H}(U_{E}^{(j)},U_{I}^{(j)}).$

  * •

Update the known values for the implicit part $U_{*}^{(i)}$, where

(3.24) $U_{*}^{(i)}=U^{n}+\Delta
t\sum^{i-1}_{j=1}a_{ij}\mathcal{H}(U_{E}^{(j)},U_{I}^{(j)}),$

and then solve

(3.25) $U_{I}^{(i)}=U_{*}^{(i)}+\Delta
ta_{ii}\mathcal{H}(U_{E}^{(i)},U_{I}^{(i)}).$

  * •

Finally, the solution $U^{n+1}$ at time level $t^{n+1}$ is accumulated by

(3.26) $U^{n+1}=U^{n}+\Delta
t\sum_{i=1}^{s}b_{i}\mathcal{H}(U_{E}^{(i)},U_{I}^{(i)}).$

In components, the procedures corresponding to $U^{(i)}_{E}$ and $U^{(i)}_{*}$
are

(3.27a) $h_{E}^{(i)}=h^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\nabla\cdot(h{\bf u})_{I}^{(j)},$ (3.27b)
$(h{\bf u})_{E}^{(i)}=(h{\bf u})^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\left(\nabla\cdot\left(\frac{h{\bf u}\otimes
h{\bf u}}{h}\right)_{E}^{(j)}+h^{(j)}_{E}\nabla H_{I,2}^{(j)}\right),$ (3.28a)
$h_{*}^{(i)}=h^{n}-\Delta t\sum_{j=1}^{i-1}a_{ij}\nabla\cdot(h{\bf
u})_{I}^{(j)},$ (3.28b) $(h{\bf u})_{*}^{(i)}=(h{\bf u})^{n}-\Delta
t\sum_{j=1}^{i-1}a_{ij}\left(\nabla\cdot\left(\frac{h{\bf u}\otimes h{\bf
u}}{h}\right)_{E}^{(j)}+h^{(j)}_{E}\nabla H_{I,2}^{(j)}\right),$

and for $U^{(i)}_{I}$ it takes the form

(3.29a) $h_{I}^{(i)}=h_{*}^{(i)}-{a}_{ii}\Delta t\nabla\cdot(h{\bf
u})_{I}^{(i)},$ (3.29b) $(h{\bf u})_{I}^{(i)}=(h{\bf
u})_{*}^{(i)}-{a}_{ii}\Delta t\left(\nabla\cdot\left(\frac{h{\bf u}\otimes
h{\bf u}}{h}\right)_{E}^{(i)}+h^{(i)}_{E}\nabla H_{I,2}^{(i)}\right).$

In order to solve the implicit components in (3.29), a similar fashion as in
the first order case can be followed. By substituting $(h{\bf u})_{I}^{(i)}$
from the second equation into the first equation, replacing $h^{(i)}_{I}$ by
$h^{(i)}_{I}=\bar{H}^{(i)}_{E}+\varepsilon^{2}H^{(i)}_{I,2}$, where
$\bar{H}^{(i)}_{E}$ is the spatial average of $h^{(i)}_{E}+b$, we obtain

(3.30) $\varepsilon^{2}H^{(i)}_{I,2}-(a_{ii}\Delta
t)^{2}\nabla\cdot\Big{(}h^{(i)}_{E}\nabla H^{(i)}_{I,2}\Big{)}=h^{**},$

with

(3.31) $h^{**}=h^{(i)}_{*}+b-\bar{H}^{(i)}_{E}-a_{ii}\Delta
t\left(\nabla\cdot(h{\bf u})^{(i)}_{*}-a_{ii}\Delta
t\nabla^{2}:\left(\frac{h{\bf u}\otimes h{\bf u}}{h}\right)^{(i)}_{E}\right).$

Lastly, the equations (3.26) can be rewritten as

(3.32a) $h^{n+1}=h^{n}-\Delta t\sum_{i=1}^{s}b_{i}\nabla\cdot(h{\bf
u})_{I}^{(i)},$ (3.32b) $(h{\bf u})^{n+1}=(h{\bf u})^{n}-\Delta
t\sum_{i=1}^{s}b_{i}\left(\nabla\cdot\left(\frac{h{\bf u}\otimes h{\bf
u}}{h}\right)_{E}^{(i)}+h^{(i)}_{E}\nabla H_{I,2}^{(i)}\right).$

For high order in space, we will adopt the finite difference WENO
reconstruction [26, 48, 49] for the first order convection terms, and central
difference for the second order and mixed derivatives. The numerical fluxes
for the convection terms are chosen in the same spirit as the first order case
described above.

We take $(hu)_{x}$ in the convection term $\nabla\cdot(h{\bf u})$ as an
example, and omit the indexes for brevity. For high order finite difference
reconstruction, the flux needs to split into an upwind and another downwind
part, e.g., for the Lax-Friedrichs flux splitting, we have

(3.33)
$(hu)^{\pm}_{i+\ell,j}=\frac{1}{2}\Big{(}(hu)_{i+\ell,j}\pm\alpha^{x}_{i,j}H_{i+\ell,j}\Big{)},\quad\ell=-r,\cdots,r,$

where $\alpha^{x}_{i,j}=\max_{h,hu}(|u|+\min(1,1/\varepsilon)\sqrt{h})$ is the
local numerical viscosity coefficient over the stencil
$S=\\{(i-r,j),\cdots,(i+r,j)\\}$. It is also important to take $H$ instead of
$h$ in (3.33), in order to preserve the still water equilibrium. The split
fluxes can be used to reconstruct $(\widehat{hu})^{\mp}_{i\pm\frac{1}{2},j}$
based on upwind and downwind WENO reconstructions. In our numerical section, a
fifth order finite difference WENO reconstruction with $r=2$ is use. The
numerical flux for $(hu)_{x}$ is defined as

(3.34)
$(\widehat{hu})_{i+\frac{1}{2},j}=(\widehat{hu})^{-}_{i+\frac{1}{2},j}+(\widehat{hu})^{+}_{i+\frac{1}{2},j}.$

The numerical flux $(\widehat{hv})_{i,j+\frac{1}{2}}$ along the $y$ direction
can be defined similarly. With these numerical fluxes, the convection term
$\nabla\cdot(h{\bf u})$ can be approximated by $\nabla_{LF}\cdot(h{\bf u})$ as
defined in (3.11). The term $\nabla\cdot\left(h{\bf u}\otimes h{\bf
u}/h\right)$ in the momentum equation can be approximated in a similar way by
the high order finite difference WENO reconstruction as $\nabla\cdot(h{\bf
u})$, e.g., for $(hu^{2})_{x}$ and $(huv)_{y}$ in the momentum equation of
$hu$, a Lax-Friedrichs flux splitting is taken as

(3.35a)
$(hu^{2})^{\pm}_{i+\ell,j}=\frac{1}{2}\Big{(}(hu^{2})_{i+\ell,j}\pm\alpha^{x}_{i,j}(hu)_{i+\ell,j}\Big{)},\quad\ell=-r,\cdots,r,$
(3.35b)
$(huv)^{\pm}_{i,j+\ell}=\frac{1}{2}\Big{(}(huv)_{i,j+\ell}\pm\alpha^{y}_{i,j}(hu)_{i,j+\ell}\Big{)},\quad\ell=-r,\cdots,r.$

Similarly for $(huv)_{x}$ and $(hv^{2})_{y}$ in the momentum equation of $hv$,
so that we get the approximation for $\nabla\cdot\left(h{\bf u}\otimes h{\bf
u}/h\right)$, which is still denoted as $\nabla_{LF}\cdot\left(h{\bf u}\otimes
h{\bf u}/h\right)$.

For the second order derivative terms appeared in $\nabla^{2}:(\frac{h{\bf
u}\otimes h{\bf u}}{h})$, a high order central difference discretization is
used, which is denoted as $\nabla^{2}_{C}:(\frac{h{\bf u}\otimes h{\bf
u}}{h})$. In our numerical section, we take a fourth order central difference
discretization. For example, along the $x$ direction, we approximate $q_{xx}$
by

$q_{xx}|_{x=x_{i}}=\frac{-q_{i-2}+16q_{i-1}-30q_{i}+16q_{i+1}-q_{i+2}}{12\Delta
x^{2}}+\mathcal{O}(\Delta x^{4}).$

For the mixed derivative term $q_{xy}$, it is discretized dimension-by-
dimension with a fourth order central difference scheme along each direction,
e.g., along the $x$ direction

$q_{x}|_{x=x_{i}}=\frac{q_{i-2}-8q_{i-1}+8q_{i+1}-q_{i+2}}{12\Delta
x^{2}}+\mathcal{O}(\Delta x^{4}).$

For the variable coefficient diffusion term $\nabla\cdot(h\nabla H)$, we take
a compact fourth order central difference discretization as developed in [5],
which is denoted as $\nabla_{C^{2}}\cdot(h\nabla H)$. Taking
$(a(x,y)q_{x})_{x}$ at the grid point $(x_{i},y_{j})$ as an example, it is
approximated by

$(a(x)q_{x})_{x}|_{(x_{i},y_{j})}=\frac{1}{\Delta x^{2}}{\bf
a}_{i,j}\begin{pmatrix}-25/144&1/3&-1/4&1/9&-1/48\\\ 1/6&5/9&-1&1/3&-1/18\\\
0&0&0&0&0\\\ -1/18&1/3&-1&5/9&1/6\\\
-1/48&1/9&-1/4&1/3&-25/144\end{pmatrix}{\bf q}^{T}_{i,j}+\mathcal{O}(\Delta
x^{4}),$

with the two vectors being

${\bf a}_{i,j}=(a_{i-2,j},a_{i-1,j},a_{i,j},a_{i+1,j},a_{i+2,j}),\quad{\bf
q}_{i,j}=(q_{i-2,j},q_{i-1,j},q_{i,j},q_{i+1,j},q_{i+2,j}).$

For high order spatial discretization, as in (3.13), we rewrite the term
$h\nabla H_{2}$ and discretize it as

(3.36)
$\nabla_{W}\big{(}\bar{H}_{E}H_{I,2}+\frac{1}{2}\varepsilon^{2}(H_{I,2})^{2}-H_{I,2}b\big{)}+H_{I,2}\nabla_{W}b,$

here $\nabla_{W}$ in the first term is a high order finite difference WENO
reconstruction, but with zero viscosity as studied in [5]. We use the same
$\nabla_{W}$ with exactly the same nonlinear weights to evaluate
$\nabla_{W}b$, for the purpose of preserving the exact still water
equilibrium. We refer to [62] for more detailed discussion of this matter.

With the above space and time discretizations, we now summarize our high order
semi-implicit scheme as follows:

  * •

First for the stage values from $i=1,\cdots,s$:

    1. (1)

update $h_{E}^{(i)}$ and $(h{\bf u})_{E}^{(i)}$ from

(3.37a) $h_{E}^{(i)}=h^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\nabla_{LF}\cdot(h{\bf u})_{I}^{(j)},$ (3.37b)
$(h{\bf u})_{E}^{(i)}=(h{\bf u})^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\left(\nabla_{LF}\cdot\left(\frac{h{\bf
u}\otimes h{\bf
u}}{h}\right)_{E}^{(j)}+{\nabla_{W}\Big{(}\bar{H}^{(j)}_{E}H^{(j)}_{I,2}+\frac{1}{2}\varepsilon^{2}(H^{(j)}_{I,2})^{2}-H^{(j)}_{I,2}b\Big{)}+H^{(j)}_{I,2}\nabla_{W}b}\right).$

    2. (2)

precompute the known values of $h_{*}^{(i)}$ and $(h{\bf u})_{*}^{(i)}$

(3.38a) $h_{*}^{(i)}=h^{n}-\Delta
t\sum_{j=1}^{i-1}a_{ij}\nabla_{LF}\cdot(h{\bf u})_{I}^{(j)},$ (3.38b) $(h{\bf
u})_{*}^{(i)}=(h{\bf u})^{n}-\Delta
t\sum_{j=1}^{i-1}a_{ij}\left(\nabla_{LF}\cdot\left(\frac{h{\bf u}\otimes h{\bf
u}}{h}\right)_{E}^{(j)}+{\nabla_{W}\Big{(}\bar{H}^{(j)}_{E}H^{(j)}_{I,2}+\frac{1}{2}\varepsilon^{2}(H^{(j)}_{I,2})^{2}-H^{(j)}_{I,2}b\Big{)}+H^{(j)}_{I,2}\nabla_{W}b}\right).$

    3. (3)

solve the linear elliptic equation to obtain $H^{(i)}_{I,2}$

(3.39a) $\varepsilon^{2}H^{(i)}_{I,2}-(a_{ii}\Delta
t)^{2}\nabla_{C^{2}}\cdot\Big{(}h^{(i)}_{E}\nabla
H^{(i)}_{I,2}\Big{)}=h^{**},$ (3.39b)
$h^{**}=h^{(i)}_{*}+b-\bar{H}^{(i)}_{E}-a_{ii}\Delta
t\left(\nabla_{LF}\cdot(h{\bf u})^{(i)}_{*}-a_{ii}\Delta
t\nabla^{2}_{C}:\left(\frac{h{\bf u}\otimes h{\bf
u}}{h}\right)^{(i)}_{E}\right).$

    4. (4)

update $h_{I}^{(i)}$ and $(h{\bf u})_{I}^{(i)}$ from

(3.40a) $h_{I}^{(i)}=h_{*}^{(i)}-{a}_{ii}\Delta t\nabla_{LF}\cdot(h{\bf
u})_{I}^{(i)},$ (3.40b) $(h{\bf u})_{I}^{(i)}=(h{\bf
u})_{*}^{(i)}-{a}_{ii}\Delta t\left(\nabla_{LF}\cdot\left(\frac{h{\bf
u}\otimes h{\bf
u}}{h}\right)_{E}^{(i)}+{\nabla_{W}\Big{(}\bar{H}^{(i)}_{E}H^{(i)}_{I,2}+\frac{1}{2}\varepsilon^{2}(H^{(i)}_{I,2})^{2}-H^{(i)}_{I,2}b\Big{)}+H^{(i)}_{I,2}\nabla_{W}b}\right).$

  * •

Update the solution at the time level $t^{n+1}$:

(3.41a) $h^{n+1}=h^{n}-\Delta t\sum_{i=1}^{s}b_{i}\nabla_{LF}\cdot(h{\bf
u})_{I}^{(i)},$ (3.41b) $(h{\bf u})^{n+1}=(h{\bf u})^{n}-\Delta
t\sum_{i=1}^{s}b_{i}\left(\nabla_{LF}\cdot\left(\frac{h{\bf u}\otimes h{\bf
u}}{h}\right)_{E}^{(i)}+{\nabla_{W}\Big{(}\bar{H}^{(i)}_{E}H^{(i)}_{I,2}+\frac{1}{2}\varepsilon^{2}(H^{(i)}_{I,2})^{2}-H^{(i)}_{I,2}b\Big{)}+H^{(i)}_{I,2}\nabla_{W}b}\right).$

###  3.3. Well-balanced property for high order semi-implicit scheme

Here we show that high order semi-implicit scheme (3.37)-(3.41) can maintain
the well-balanced property for the still water equilibrium (1.2). We have the
following theorem:

######  Theorem 3.1.

The high order semi-implicit scheme (3.37)-(3.41) is well-balanced for the
still water equilibrium, in the sence that, if initially the water is at
still, namely $H^{0}=$ Const. and $(h{\bf u})^{0}={\bf 0}$, the scheme can
maintain still water at any later time with $H^{n}=$ Const. and $(h{\bf
u})^{n}={\bf 0}$.

###### Proof.

We prove this theorem using the mathematical induction. Assume at the time
step $t^{n}$, we have $H^{n}=$ Const. and ${\bf u}^{n}={\bf 0}$. First for
$U^{(0)}_{E}=U^{(0)}_{I}=U^{n}$, similar to the discussion in the first order
semi-implicit scheme, we have

(3.42) $\nabla_{LF}\cdot(h{\bf
u})_{I}^{n}=0,\quad\nabla_{LF}\cdot\left(\frac{h{\bf u}\otimes h{\bf
u}}{h}\right)_{E}^{n}=0,\quad{H^{n}_{I,2}=0},\quad\nabla^{2}_{C}:\left(\frac{h{\bf
u}\otimes h{\bf u}}{h}\right)^{n}=0.$

By the induction hypothesis, we assume $H^{(j)}_{E}=H^{(j)}_{I}=$ Const. and
$(h{\bf u})^{(j)}_{E}=(h{\bf u})^{(j)}_{I}={\bf 0}$ hold for any $j\leq i-1$,
from which we have for $j=1,\cdots,i-1$:

(3.43) $\nabla_{LF}\cdot(h{\bf
u})_{I}^{(j)}=0,\quad\nabla_{LF}\cdot\left(\frac{h{\bf u}\otimes h{\bf
u}}{h}\right)_{E}^{(j)}=0,\quad{H^{(j)}_{I,2}=0},\quad\nabla^{2}_{C}:\left(\frac{h{\bf
u}\otimes h{\bf u}}{h}\right)^{(j)}_{E}=0.$

The goal is to show that $H^{(i)}_{E}=H^{(i)}_{I}=$ Const. and $(h{\bf
u})^{(i)}_{E}=(h{\bf u})^{(i)}_{I}={\bf 0}$. From (3.37), we have
$h^{(i)}_{E}=h^{n}$ so that $H^{(i)}_{E}=h^{n}+b=$ Const., and $(h{\bf
u})^{(i)}_{E}=(h{\bf u})^{n}={\bf 0}$. Similarly $h_{*}^{(i)}+b=h^{n}+b=$
Const., and $(h{\bf u})^{(i)}_{*}=(h{\bf u})^{n}={\bf 0}$ from (3.38). With
these, we conclude that $h^{**}=0$ from (3.39b), and solving the elliptic
equation (3.39a) with a positive definite matrix leads to $H^{(i)}_{I,2}=0$.
Furthermore, due to $(h{\bf u})^{(i)}_{E}={\bf 0}$ and $H^{(i)}_{I,2}=0$, we
get $(h{\bf u})^{(i)}_{I}=(h{\bf u})^{(i)}_{*}={\bf 0}$ from (3.40b). It
follows that $h^{(i)}_{I}=h^{(i)}_{*}$ from (3.40a), so that
$H^{(i)}_{I}=h^{(i)}_{*}+b=$ Const., and we complete the mathematical
induction.

Since (3.43) holds for $j=1,\cdots,s$, substituting them into (3.41), we
obtain $H^{n+1}=H^{n}=$ Const. and $(h{\bf u})^{n+1}=(h{\bf u})^{n}={\bf 0}$.
Therefore, the well-balanced property is preserved and this finishes the
proof. ∎

###  3.4. Asymptotic preserving and asymptotically accurate properties

In this section, we formally prove the AP property for the first order semi-
implicit scheme (3.8), and the AA property for the high order semi-implicit
scheme (3.27)-(3.32). When we discuss the AP or AA property, we focus on the
time discretization while keeping the space continuous. First we have the
following theorem.

######  Theorem 3.2.

The first order semi-implicit scheme (3.8) with space continuous is asymptotic
preserving, in the sense that, with no-slip or periodic boundary condition, at
the leading order asymptotic expansions, the scheme (3.8) is a consistent
approximation of the lake equations (2.10) at the zero Froude number limit.

###### Proof.

To prove the theorem, we assume the following expansions of the solutions at
all time levels, i.e., $h^{n}({\bf x}):=h({\bf x},t^{n})$ and ${\bf
u}^{n}({\bf x}):={\bf u}({\bf x},t^{n})$ admit

(3.44) $h^{n}({\bf x})=h^{n}_{0}({\bf x})+\varepsilon^{2}H^{n}_{2}({\bf
x}),\quad{\bf u}^{n}({\bf x})={\bf u}^{n}_{0}({\bf x})+\varepsilon{\bf
u}^{n}_{1}({\bf x}),$

and correspondingly the water surface level $H^{n}({\bf x}):=H({\bf
x},t^{n})=h({\bf x},t^{n})+b({\bf x})$ takes the form

(3.45) $H^{n}({\bf x})=h^{n}_{0}({\bf x})+b({\bf
x})+\varepsilon^{2}H^{n}_{2}({\bf x}),$

where $H_{0}=h^{n}_{0}({\bf x})+b({\bf x})=$ Const., namely $h^{n}_{0}({\bf
x})$ does not depend on $n$.

We plug them into the semi-discrete scheme (3.8), with the first equation
equivalent to the first equation of (3.1). From
$H^{n+1}=\bar{H}^{n}+\varepsilon^{2}H^{n+1}_{2}$, it yields
$\bar{H}^{n}=H_{0}=$ Const. Equating to zero for the
$\mathcal{O}(\varepsilon^{0})$ terms, we have

(3.46) $\left\\{\begin{aligned} &\nabla\cdot(h_{0}{\bf u}_{0})^{n+1}=0,\\\
&\frac{(h_{0}{\bf u}_{0})^{n+1}-(h_{0}{\bf u}_{0})^{n}}{\Delta
t}+\nabla\cdot\left(\frac{h_{0}{\bf u}_{0}\otimes h_{0}{\bf
u}_{0}}{h_{0}}\right)^{n}+h_{0}^{n}\nabla H^{n+1}_{2}=0,\end{aligned}\right.$

which is a consistent discretization to the lake equations (2.10), with
$H^{n+1}_{2}$ solved from

(3.47) $-\nabla\cdot(h_{0}^{n}\nabla
H^{n+1}_{2})=\nabla^{2}:\left(\frac{h_{0}{\bf u}_{0}\otimes h_{0}{\bf
u}_{0}}{h_{0}}\right)^{n}.$

∎

Now we are ready to present the AA property for the high order semi-implicit
scheme (3.27)-(3.32), that is, the scheme maintains its temporal order of
accuracy for the lake equations (2.10) at the zero Froude limit when
$\varepsilon\rightarrow 0$ [45, 6]. To have the AA property, it is crucial
that the IMEX Runge-Kutta scheme (3.21) is stiffly accurate (SA), namely the
implicit part satisfies $b^{T}={\bf e}^{T}_{s}A$, where ${\bf
e}^{T}_{s}=(0,\cdots,0,1)$ [4, 6]. Besides, the initial conditions
$(h^{0}({\bf x}),h^{0}({\bf x}){\bf u}^{0}({\bf x}))$ need to be well-
prepared, in the sense that

(3.48) $h^{0}({\bf x})=h_{0}({\bf x})+\varepsilon^{2}H^{0}_{2}({\bf
x}),\quad{\bf u}^{0}({\bf x})={\bf u}_{0}({\bf x})+\varepsilon{\bf u}_{1}({\bf
x}),\quad\text{ and }\quad h_{0}({\bf x})+b({\bf x})=H_{0}=\text{Const}.$

We have the following theorem about the AA property of the high order semi-
implicit scheme.

######  Theorem 3.3.

For the high order semi-implicit scheme (3.27)-(3.32) of temporal order $p$,
when applied to the system (2.1) on a bounded domain with no-slip or periodic
boundary condition, suppose the IMEX Runge-Kutta scheme (3.21) is stiffly
accurate, and the initial conditions $(h^{0}({\bf x}),h^{0}({\bf x}){\bf
u}^{0}({\bf x}))$ are well prepared (3.48). Denoting by ${\bf V}^{1}({\bf
x};\varepsilon)=(h^{1}({\bf x};\varepsilon),h^{1}({\bf x};\varepsilon){\bf
u}^{1}({\bf x};\varepsilon))$ the numerical solution after one time step, we
have

(3.49) $\lim_{\varepsilon\to 0}h^{1}({\bf x};\varepsilon)+b({\bf
x})=H_{0},\quad\lim_{\varepsilon\to 0}\nabla\cdot\Big{(}h^{1}({\bf
x};\varepsilon){\bf u}^{1}({{\bf x}};\varepsilon)\Big{)}=0.$

Furthermore, let ${\bf V}_{lake}({\bf x},t)=(h_{lake}({\bf x},t),h_{lake}({\bf
x},t){\bf u}_{lake}({\bf x},t))$ be the exact solution of the lake equations
(2.10) with initial conditions $(h^{0}({\bf x}),h^{0}({\bf x}){\bf u}^{0}({\bf
x}))$, one has the one-step error estimate

(3.50) $\lim_{\varepsilon\to 0}{\bf V}^{1}({\bf x};\varepsilon)={\bf
V}_{lake}({\bf x},\Delta t)+\mathcal{O}(\Delta t^{p+1}),$

i.e., the high order semi-implicit scheme is AA.

The proof follows from the same structure as in [6] by the mathematical
induction, and is skipped here.

##  4\. Numerical tests

In this section, we will perform some numerical tests with the Froude number
ranging from $0$ to ${\mathcal{O}}(1)$. The fifth order finite difference WENO
reconstruction [48, 49, 62] is used for the first order spatial derivatives,
and the fourth order (compact) central difference discretizations for the
second order derivatives. In time we employ a third order SA IMEX Runge-Kutta
scheme SI-IMEX(4,4,3) from [6], with the double Butcher tableau given by

Explicit : (4.6) $\displaystyle\begin{array}[]{c|cccc}0&0&0&0&0\\\
\gamma&\gamma&0&0&0\\\ 0.717933260754&1.243893189483&-0.525959928729&0&0\\\
1&0.630412558153&0.786580740199&-0.416993298352&0\\\ \hline\cr
0&0&1.208496649176&-0.644363170684&\gamma\end{array},$ (4.7) Implicit : (4.13)
$\displaystyle\begin{array}[]{c|cccc}\gamma&\gamma&0&0&0\\\
\gamma&0&\gamma&0&0\\\ 0.717933260754&0&0.282066739245&\gamma&0\\\
1&0&1.208496649176&-0.644363170684&\gamma\\\
\hline\cr&0&1.208496649176&-0.644363170684&\gamma\end{array},$

where $\gamma=0.435866521508$.

The time step is taken as

$\Delta t=\text{CFL}\,\Delta x/\Lambda,\qquad\Lambda=\max\\{|{\bf
u}|+\min(1,1/\varepsilon)\sqrt{h}\\},$

and $\text{CFL}=0.2$ is used. For the accuracy tests, the time step is
modified to $\Delta t=\text{CFL}\Delta x^{5/3}/\Lambda$, for better
observation of spatial orders. $N$ or $N^{2}$ uniform gird points are used for
1D and 2D problems respectively, except otherwise specified.

For the Froude number $\varepsilon$ of ${\mathcal{O}}(1)$, e.g.
$\varepsilon=\frac{1}{\sqrt{g}}$ and $g=9.812$ is the gravitational constant,
we will compare our results to reference solutions, which are produced by the
fifth order well-balanced finite difference WENO scheme developed by Xing and
Shu [62]. We refer it as “WB-Xing” in the following.

###  4.1. One dimensional case

######  Example 4.1.

(Accuracy test) We first consider an example with smooth initial conditions
and a non-zero bottom topology, which are given by [62]

(4.14) $h(x,0)=5+\exp(\cos(2\pi x)),\quad(hu)(x,0)=\sin(\cos(2\pi x)),\quad
b(x)=\sin^{2}(\pi x),\quad x\in[0,1].$

Periodic boundary condition is used, with the Froude number
$\varepsilon=\frac{1}{\sqrt{g}}$. We take mesh grid points $N=40\times 2^{i}$
for $i=1,\cdots,5$. Since the exact solution is not available, numerical
errors are computed by comparing the numerical solutions at two successive
mesh grid points from refinement. For this case, the errors are computed for
the momentum $hu$, at a final time $T=0.1$. The results are shown in Table
4.1, we can see that it approaches fifth order, which is similar to the
results in [62].

Table 4.1. Example 4.1: numerical errors and orders of accuracy for the momentum $hu$. $T=0.1$. N | 80 | 160 | 320 | 640 | 1280
---|---|---|---|---|---
$L_{1}$ error | 3.35E-02 | 4.61E-03 | 4.44E-04 | 2.06E-05 | 6.99E-07
order | – | 2.86 | 3.37 | 4.43 | 4.88

######  Example 4.2.

(Accuracy test for a range of $\varepsilon$) In this example, we try to test
the orders of accuracy for our scheme in different regimes of the Froude
number. We take the non-flat bottom topological function $b(x)$ as

$b(x)=1+\sin(2\pi x),$

with initial conditions

(4.15) $\displaystyle h(x,0)=10-b(x)+\varepsilon^{2}\exp(\sin(2\pi x)),$
$\displaystyle(hu)(x,0)=1+\varepsilon^{2}\sin(2\pi x).$

Similarly, periodic boundary condition is used. Numerical errors are computed
in the same way as in the previous example. Three different Froude numbers
$\varepsilon=1,10^{-2},10^{-4}$ are taken, with a final time $T=0.05$. And the
computational domain is $x\in[0,2]$. Numerical errors and orders are shown in
Table 4.2. From this table, we can see that for all $\varepsilon$’s, almost
fifth order accuracy can be observed, which indicate that our scheme is
asymptotically accurate.

Table 4.2.  Example 4.2: numerical errors and orders of accuracy for the momentum $hu$ with initial conditions (4.15). $T=0.05$. $\varepsilon$ | N | 80 | 160 | 320 | 640
---|---|---|---|---|---
1 |  $L^{1}$ error | 6.09E-03 | 3.23E-04 | 1.16E-05 | 4.05E-07
order | – | 4.24 | 4.80 | 4.84
$10^{-2}$ |  $L^{1}$ error | 8.27E-03 | 3.75E-04 | 2.70E-05 | 1.06E-06
order | – | 4.46 | 3.79 | 4.67
$10^{-4}$ |  $L^{1}$ error | 4.58E-05 | 4.92E-06 | 1.26E-06 | 5.58E-08
order | – | 3.22 | 1.97 | 4.49

######  Example 4.3.

(A small perturbation of a steady-state water) This example was first proposed
by LeVeque in [35] and later studied by Xing and Shu in [62]. There is a small
perturbation on a quasi-stationary water, moving over a non-flat bottom
topography. The bottom function is smooth, which is given by

(4.16) $b(x)=\left\\{\begin{aligned} &0.25(\cos(10\pi(x-1.5))+1),&\text{if
}1.4\leq x\leq 1.6;\\\ &0,&\text{otherwise},\end{aligned}\right.$

and the initial conditions are:

(4.17a) $h(x,0)=\left\\{\begin{aligned} &1-b(x)+\eta,&\text{if }1.1\leq x\leq
1.2;\\\ &1-b(x),&\text{otherwise},\end{aligned}\right.$ (4.17b) $hu(x,0)=0,$

on the domain $x\in[0,2]$, see Fig. 4.1. $\eta$ is the magnitude of
perturbation. Two cases are considered: $\eta=0.2$ (big pulse) and
$\eta=0.001$ (small pulse). The Froude number is taken to be
$\varepsilon=\frac{1}{\sqrt{g}}$. After the perturbation moves over the non-
flat bottom, two disturbances will generate and one propagates to the left and
the other to the right, both with a speed $\sqrt{gh}$.

This example is used to test the well-balanced property of the numerical
scheme. For non well-balanced schemes, numerical errors may pollute the small
perturbations. We show the water surface level $H$ and momentum $hu$, for
$\eta=0.2$ and $\eta=0.001$ in Fig. 4.2 and Fig. 4.3 respectively, at a final
time $T=0.2$ with $N=200$. We compare the solutions to the reference solutions
of “WB-Xing” with $N=3000$. It can be observed that, for both cases, our
solutions can well capture the disturbances and match the reference solutions.

(a) $\eta=0.2$

(b) $\eta=0.001$

Figure 4.1.  The initial water surface level $h+b$ (4.17a) and the bottom $b$
(4.16) for Example 4.3. Left: $\eta=0.2$; Right:$\eta=0.001$.

Figure 4.2. The water surface level $h+b$ (left) and momentum $hu$ (right) at
time $T=0.2$ with $\eta=0.2$ for Example 4.3.

Figure 4.3. The water surface level $h+b$ (left) and momentum $hu$ (right) at
time $T=0.2$ with $\eta=0.001$ for Example 4.3.

######  Example 4.4.

(Dam breaking) The dam breaking problem over a rectangular bump is widely used
to test the oscillation-free property of numerical schemes for the shallow
water equations [58]. The bottom function is defined as

(4.18) $b(x)=\left\\{\begin{aligned} &8,&\text{if }|x-750|\leq 1500/8;\\\
&0,&\text{otherwise},\end{aligned}\right.$

and the initial conditions are:

(4.19) $(hu)(x,0)=0\quad\text{ and }\quad h(x,0)=\left\\{\begin{aligned}
&20-b(x),&\text{if }x\leq 750;\\\
&15-b(x),&\text{otherwise},\end{aligned}\right.$

on a computational domain $x\in[0,1500]$. The Froude number is taken as
$\varepsilon=\frac{1}{\sqrt{g}}$. The inflow and outflow boundary conditions
are set the same as the initial values on the left and right respectively. As
time evolves, the initial jump on $h$ will generate two waves. One is a
rarefaction wave traveling to the left, and the other is a shock traveling to
the right. For this example, we show the water surface level $h+b$ on the mesh
points $N=500$ at two different times $T=15$ and $60$ in Fig. 4.4 and Fig.
4.5, respectively. We also compare them to the reference solutions on $N=3000$
with “WB-Xing” method. The results match each other well.

Figure 4.4. The water surface level for the dam breaking problem of Example
4.4 at $T=15$. Left: initial surface level $h+b$, final surface level $h+b$
and bottom topology $b$; Right: comparisons to the reference solutions.

Figure 4.5. The water surface level for the dam breaking problem of Example
4.4 at $T=60$. Left: initial surface level $h+b$, final surface level $h+b$
and bottom topology $b$; Right: comparison to the reference solutions.

######  Example 4.5.

(Lake at rest) In this example, we consider a still water initially to test
the well-balanced property of our scheme. We take a non-smooth bottom topology
given by

(4.20) $b(x)=\left\\{\begin{aligned} &4,&\text{if }4\leq x\leq 8;\\\
&0,&\text{otherwise},\end{aligned}\right.$

and the initial conditions are

(4.21) $(h+b)(x,0)=10,\quad(hu)(x,0)=0,$

on a computational domain $[0,10]$ with the Froude number set to be
$\varepsilon=\frac{1}{\sqrt{g}}$. Periodic boundary condition is adopted.

In Fig. 4.6, we show the variation of the water surface level, and the
momentum at final time $T=10$. We can clearly see the errors are within
machine precision, namely, the still water equilibrium is well preserved.

To further show the ability of our scheme, we set the initial velocity as
$u(x,0)=1$ for a moving water. Due to the non-flat bottom $b$, the initial
water equilibrium will be destroyed. We compute the numerical solution to
$T=0.1$. The results are shown in Fig. 4.7 and compared to the reference
solutions from “WB-Xing”. We can see the results still match each other, and
the discontinuities are well captured without any artificial oscillation.

Figure 4.6. The lake at rest with non-smooth bottom for Example 4.5 at $T=10$.
Left: the disturbance of $H$; Right: the numerical result for momentum $hu$.

Figure 4.7. Moving water with non-smooth bottom for Example 4.5 at $T=0.1$.
Left: the water surface level $h+b$; Right: the momentum $hu$.

######  Example 4.6.

In this example, we would like to test our scheme for a multiscale wave, which
was studied in [30, 37]. The initial conditions are taken as

(4.22a) $H(x,0)=1+\frac{\sigma(x)}{2}\sin(\varepsilon 40\pi
x)+\varepsilon(1+\cos(\varepsilon\pi x)),$ (4.22b)
$u(x,0)=\sqrt{2}(1+\cos(\varepsilon\pi x)),\quad\quad b(x)=0,$

with the Froude number $\varepsilon=0.02$, and

$\sigma(x)=\left\\{\begin{aligned} &0.5(1-\cos(0.1\pi x)),&\text{if }0\leq
x\leq 20;\\\ &0,&\text{otherwise}.\end{aligned}\right.$

The computational domain is $[-51,51]$ with periodic boundary condition.

We show the numerical solutions of the multiscale wave propagation in Fig. 4.8
on a uniform mesh of $N=2040$. We compare our results to those produced by the
explicit scheme of “WB-Xing” on the same mesh. We can see that under this mesh
size, the results of our AP scheme match those from the explicit “WB-Xing”
scheme very well.

(a) $t=0$

(b) $t=0.2$

(c) $t=0.5$

(d) $t=1.0$

(e) $t=2.4$

(f) $t=4.1$

Figure 4.8.  The numerical solution of water surface level $h+b$ for Example
4.6 with a uniform mesh of $N=2040$.

###  4.2. Two-dimensional case

######  Example 4.7.

(Accuracy test) For this 2D example, we consider a smooth non-flat bottom
function to be

(4.23) $b(x,y)=\sin(2\pi x)+\cos(2\pi y)+2,$

and the initial conditions are

(4.24) $\left\\{\begin{aligned} h(x,y,0)&=10-b(x,y)+\varepsilon^{2}\sin(2\pi
x)\cos(2\pi y),\\\ (hu)(x,y,0)&=\sin(2\pi x)\cos(2\pi y),\\\
(hv)(x,y,0)&=-\cos(2\pi x)\sin(2\pi y),\end{aligned}\right.$

on a computational domain $[0,1]^{2}$ with periodic boundary conditions in
both directions. Note that the initial conditions (4.24) are set to be well-
prepared (3.48).

We take three different $\varepsilon$’s of $\varepsilon=1,10^{-2}$ and
$10^{-4}$. We compute the solution up to a final time $T=0.05$ on mesh grid
points of $N^{2}$. Since the exact solution is not available, the numerical
errors are computed by comparing the solutions on two consecutive meshes. The
$L^{1}$ errors and orders of accuracy are shown in Table 4.3. We can see high
order accuracy can be obtained for all Froude numbers.

Table 4.3.  Example 4.7. Numerical errors and orders of accuracy for the momentum $hu$ with Froude number $\varepsilon=1,10^{-2},10^{-4}$. $\varepsilon$ | N | 16 | 32 | 64 | 128 | 256
---|---|---|---|---|---|---
1 |  $L^{1}$ error | 6.56E-02 | 3.74E-03 | 1.36E-04 | 4.39E-06 | 1.39E-07
order | – | 4.13 | 4.78 | 4.96 | 4.98
$10^{-2}$ |  $L^{1}$ error | 2.47E-02 | 1.46E-03 | 1.69E-04 | 6.11E-06 | 1.91E-07
order | – | 4.08 | 3.11 | 4.79 | 5.00
$10^{-4}$ |  $L^{1}$ error | 2.63E-02 | 1.33E-03 | 4.88E-05 | 1.63E-06 | 1.09E-07
order | – | 4.31 | 4.76 | 4.90 | 3.91

######  Example 4.8.

(A small perturbation of 2D steady-state water) For this example, we try to
test our scheme for the capability of capturing the perturbation on a
stationary water in the two dimensional case, which has been studied in [35,
62].

The bottom topography is an isolated elliptical shaped hump

(4.25) $b(x,y)=0.8e^{-5(x-0.9)^{2}-50(y-0.5)^{2}}$

and the initial conditions are

(4.26a) $h(x,y,0)=\left\\{\begin{aligned} &1-b(x,y)+0.01,&\text{if}\quad
0.05\leq x\leq 0.15;\\\ &1-b(x,y),&\text{otherwise}.\end{aligned}\right.$
(4.26b) $hu=hv=0,$

on the computational domain $[0,2]\times[0,1]$, with outflow boundary in the
$x$ direction and periodic boundary in the $y$ direction. The Froude number is
set as $\varepsilon=\frac{1}{\sqrt{g}}$. We show the numerical results of
surface level $H=h+b$ on two different meshes $200\times 100$ and $400\times
200$ in Fig. 4.9. The initial perturbation is separated into two wave
propagating to the left and right. With the left-propagating wave moving out
of the domain, the right-propagating wave interacts with the non-flat bottom
topography, and is well captured by the proposed method. We can observe that
the numerical results are comparable to those of “WB-Xing” method in [62].

(a) surface level at $t=0.12$

(b) surface level at $t=0.12$

(c) surface level at $t=0.24$

(d) surface level at $t=0.24$

(e) surface level at $t=0.36$

(f) surface level at $t=0.36$

(g) surface level at $t=0.48$

(h) surface level at $t=0.48$

(i) surface level at $t=0.6$

(j) surface level at $t=0.6$

Figure 4.9.  Numerical solutions of the surface level $H=h+b$ for Example 4.8.
From top to bottom: at $t=0.12$ from 0.9993 to 1.0059; at $t=0.24$ from 0.9980
to 1.0115; at $t=0.36$ from 0.918 to 1.00872; $t=0.48$ from 0.9912 to 1.0044;
$t=0.6$ from 0.9961 to 1.00432. 30 contour lines are used. Left: $200\times
100$ uniform mesh. Right: $400\times 200$ uniform mesh.

######  Example 4.9.

(Traveling vortex) Now we consider a traveling vortex in the two dimensional
case [2]. The computational domain is $[0,2]\times[0,1]$, and the initial
conditions are given by

(4.27a) $H(x,y,0)=110+\left\\{\begin{aligned}
&\left(\frac{\varepsilon\Gamma}{\omega}\right)^{2}(k(\omega\Gamma_{c})-k(\pi)),\quad&\text{if
}\omega\Gamma_{c}\leq\pi;\\\ &0,\quad&\text{otherwise},\end{aligned}\right.$
(4.27b) $u(x,y,0)=2+\left\\{\begin{aligned}
&\Gamma(1+\cos(\omega\Gamma_{c}))(0.5-y),\quad&\text{if
}\omega\Gamma_{c}\leq\pi;\\\ &0,\quad&\text{otherwise},\end{aligned}\right.$
(4.27c) $v(x,y,0)=\left\\{\begin{aligned}
&\Gamma(1+\cos(\omega\Gamma_{c}))(x-0.5),\quad&\text{if
}\omega\Gamma_{c}\leq\pi;\\\ &0,\quad&\text{otherwise},\end{aligned}\right.$

where

(4.28)
$\Gamma_{c}=\sqrt{(x-0.5)^{2}+(y-0.5)^{2}},\quad\quad\Gamma=8,\quad\omega=4\pi,$

and

(4.29)
$k(\xi)=2\cos(\xi)+2\xi\sin(\xi)+\frac{1}{8}\cos(2\xi)+\frac{\xi}{4}\sin(2\xi)+\frac{3}{4}\xi^{2}.$

The center of the vortex is initially located at $(0.5,0.5)$, and then
propagates with a speed $u_{ref}=2$ along the horizontal direction. Periodic
boundary conditions are used. For a flat bottom, the vortex could be referred
as traveling only along the x-direction, where the exact solutions are given
as follows [47]

(4.30) $H(x,y,t)=H(u-2t,y,0),\quad u(x,y,t)=u(x-2t,y,0),\quad
v(x,y,t)=v(x-2t,y,0).$

Note that the velocity can be decomposed as ${\bf u}={\bf u}_{ref}+{\bf
u}^{\prime}$, where ${\bf u}_{ref}$ is the background traveling velocity and
${\bf u}^{\prime}$ is the rotating part which satisfies

$\nabla\cdot{\bf u}^{\prime}=0,\quad({\bf u}^{\prime}\cdot\nabla){\bf
u}^{\prime}+\nabla H=0.$

Namely, the rotating part ${\bf u}^{\prime}$ is divergence free and balanced
with $\nabla H$, so it performs as local self-rotating. We show the numerical
solutions on a mesh gird of $200\times 100$ at the final time $T=1$ in Fig.
4.10 and Fig. 4.11, and comparing our results to those produced by the
explicit scheme of “WB-Xing” on the same mesh, for three different Froude
numbers $\varepsilon=1,0.05,0.01$. The perturbation of the water surface level
$H$ from a constant level $110$ is at the scale of $\varepsilon^{2}$. We can
see that for large Froude number $\varepsilon=1$, both schemes capture the
traveling wave well. However, as the Froude number becomes small, e.g.
$\varepsilon=0.05$, our AP scheme can still keep the good shape of the vortex,
while the results from the explicit “WB-Xing” scheme have been greatly damped,
due to large numerical viscosities which are inversely proportional to the
Froude number $\varepsilon$. For the case of $\varepsilon=0.01$, our AP scheme
still has good performance, while the wave has been totally damped out for the
“WB-Xing” scheme, and numerical noises from the damped wave spreading up to
the boundary now pollute the whole computational domain, which is also the
case for smaller $\varepsilon$’s.

(a) $\varepsilon=1$

(b) $\varepsilon=1$

(c) $\varepsilon=0.05$

(d) $\varepsilon=0.05$

(e) $\varepsilon=0.01$

(f) $\varepsilon=0.01$

Figure 4.10.  Example 4.9. Numerical solutions about surface level for the
traveling vortex at time $T=1$, on a mesh $200\times 100$. $H^{\prime}=H-110$
is the deviation from the water surface level of $110$. From top to bottom
$\varepsilon=1,0.05,0.01$ respectively. Left: IMEX; Right: WB-Xing.

(a) $\varepsilon=1$

(b) $\varepsilon=1$

(c) $\varepsilon=0.05$

(d) $\varepsilon=0.05$

(e) $\varepsilon=0.01$

(f) $\varepsilon=0.01$

Figure 4.11.  Example 4.9. Numerical solutions about momentum $hu$ for the
traveling vortex at time $T=1$, on a mesh $200\times 100$. From top to bottom
$\varepsilon=1,0.05,0.01$ respectively. Left: IMEX; Right: WB-Xing.

Next we add a non-flat bottom which is variant in the $x$ direction,

$b(x,y)=e^{-5(x-1)^{2}},$

and keep others the same as in (4.27). In this case, the water surface level
would be perturbed a little due to the non-flat bottom, but the vortex still
travels almost the same. A similar example has been studied in [2, 37]. In
Fig. 4.12, we show the numerical solutions at several different times
$T=0,0.3,0.6,1.0$ with $\varepsilon$ chosen as $0.05$. The traveling vortex
can also be well captured in this case. Similarly the solutions of the “WB-
Xing” scheme have been damped.

(a) $t=0$

(b) $t=0$

(c) $t=0.3$

(d) $t=0.3$

(e) $t=0.6$

(f) $t=0.6$

(g) $t=1.0$

(h) $t=1.0$

Figure 4.12.  Example 4.9. Numerical solutions of the traveling vortex with a
nonflat bottom, on a mesh $200\times 100$. $H^{\prime}=H-110$ is the deviation
from the water surface level. From top to bottom $T=0,0.3,0.6,1.0$
respectively, and $\varepsilon=0.05$. Left: IMEX; Right: WB-Xing.

Finally, in Table 4.4 we compare the CPU cost for the two schemes with
different $\varepsilon$’s, for the cases considered above. We can find that,
the CPU time of the IMEX scheme is less than the explicit “WB-Xing” scheme,
especially in the low Froude regime, as the explicit method requires a much
smaller time step for stability. Generally, the IMEX scheme would be much more
efficient than the explicit one in the low Froude regime.

Table 4.4.  Example 4.9. The CPU time (seconds) for two schemes with different Froude numbers, with flat and non-flat bottom topographies.  Bottom topology | $\varepsilon$ | IMEX | WB-Xing
---|---|---|---
$b=0$ | 1 | 4071.7 | 5658.3
0.05 | 5342.3 | 79818.8
0.01 | 12444.2 | 378733.6
$b\neq 0$ | 0.05 | 6634.5 | 81090.6

##  5\. Conclusion

In this paper, a high order semi-implicit asymptotic preserving scheme for the
shallow water equations with a non-flat bottom topography is developed. The
scheme is shown to be well-balanced, asymptotic preserving and asymptotically
accurate. Numerical results in 1D and 2D have demonstrated the well-balanced
property, the capability of capturing small perturbations of still water
equilibrium, high order accuracy and asymptotic preserving for all ranges of
Froude numbers. As compared to the explicit “WB-Xing” scheme, the semi-
implicit AP scheme performs almost the same for large Froude numbers while
capturing small perturbations well, and is in general much more efficient in
the low Froude regime.

## References

  * [1] E. Audusse, F. Bouchut, M.-O. Bristeau, R. Klein, and B. Perthame.  A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows.  SIAM Journal on Scientific Computing, 25(6):2050–2065, 2004.
  * [2] G. Bispen, K. R. Arun, M. Lukáčová-Medvidová, and S. Noelle.  IMEX large time step finite volume methods for low Froude number shallow water flows.  Communications in Computational Physics, 16(2):307–347, 2014.
  * [3] G. Bispen, M. Lukáčová-Medvid’ová, and L. Yelash.  Asymptotic preserving IMEX finite volume schemes for low Mach number Euler equations with gravitation.  Journal of Computational Physics, 335:222–248, 2017.
  * [4] S. Boscarino, L. Pareschi, and G. Russo.  Implicit-explicit runge–kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit.  SIAM Journal on Scientific Computing, 35(1):A22–A51, 2013.
  * [5] S. Boscarino, J.-M. Qiu, G. Russo, and T. Xiong.  A high order semi-implicit IMEX WENO scheme for the all-Mach isentropic Euler system.  Journal of Computational Physics, 392:594–618, 2019.
  * [6] S. Boscarino, J.-M. Qiu, G. Russo, and T. Xiong.  High order semi-implicit WENO schemes for all Mach full Euler system of gas dynamics.  submitted.
  * [7] S. Boscarino, G. Russo, and L. Scandurra.  All Mach number second order semi-implicit scheme for the Euler equations of gas dynamics.  Journal of Scientific Computing, 77(2):850–884, 2018.
  * [8] S. Busto, L. Río-Martín, M. E. Vázquez-Cendón, and M. Dumbser.  A semi-implicit hybrid finite volume/finite element scheme for all Mach number flows on staggered unstructured meshes.  Applied Mathematics and Computation, 402:126117, 2021.
  * [9] J. C. Butcher.  Numerical Methods for Ordinary Differential Equations.  Third Edition, John Wiley & Sons Ltd, 2016.
  * [10] M. Castro, J. Gallardo, and C. Parés.  High order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. Applications to shallow-water systems.  Mathematics of computation, 75(255):1103–1134, 2006.
  * [11] S. Chen, B. Lin, Y. Li, and C. Yan.  HLLC+: Low-Mach Shock-stable HLLC-type Riemann solver for all-speed flows.  SIAM Journal on Scientific Computing, 42(4):B921–B950, 2020.
  * [12] F. Cordier, P. Degond, and A. Kumbaro.  An asymptotic-preserving all-speed scheme for the Euler and Navier-Stokes equations.  Journal of Computational Physics, 231(17):5685–5704, 2012.
  * [13] F. Couderc, A. Duran, and J.-P. Vila.  An explicit asymptotic preserving low Froude scheme for the multilayer shallow water model with density stratification.  Journal of Computational Physics, 343:235–270, 2017.
  * [14] P. Degond and M. Tang.  All speed scheme for the low mach number limit of the isentropic euler equations.  Communications in Computational Physics, 10(1):1–31, 2011.
  * [15] S. Dellacherie.  Analysis of Godunov type schemes applied to the compressibleEuler system at low Mach number.  Journal of Computational Physics, 229(4):978–1016, 2010.
  * [16] F. Denner, F. Evrard, and B. van Wachem.  Conservative finite-volume framework and pressure-based algorithm for flows of incompressible, ideal-gas and real-gas fluids at all speeds.  Journal of Computational Physics, 409:109348, 2020.
  * [17] F. Denner, C.-N. Xiao, and B. van Wachem.  Pressure-based algorithm for compressible interfacial flows with acoustically-conservative interface discretization.  Journal of Computational Physics, 367:192–234, 2018.
  * [18] G. Dimarco, R. Loubere, and M.-H. Vignal.  Study of a new asymptotic preserving scheme for the euler system in the low mach number limit.  SIAM journal on Scientific Computing, 39(5):A2099–A2128, 2017.
  * [19] A. Duran, F. Marche, R. Turpault, and C. Berthon.  Asymptotic preserving scheme for the shallow water equations with source terms on unstructured meshes.  Journal of Computational Physics, 287:184–206, 2015.
  * [20] Z. Gao and G. Hu.  High order well-balanced weighted compact nonlinear schemes for shallow water equations.  Communications in Computational Physics, 22(4):1049–1068, 2017\.
  * [21] F. X. Giraldo, J. F. Kelly, and E. M. Constantinescu.  Implicit-explicit formulations of a three-dimensional nonhydrostatic unified model of the atmosphere (numa).  SIAM Journal on Scientific Computing, 35(5):B1162–B1194, 2013.
  * [22] F. X. Giraldo and M. Restelli.  High-order semi-implicit time-integrators for a triangular discontinuous Galerkin oceanic shallow water model.  International Journal for Numerical Methods in Fluids, 63(9):1077–1102, 2010.
  * [23] H. P. Greenspan.  The theory of rotating fluids.  CUP Archive, 1968.
  * [24] J. Haack, S. Jin, and J.-G. Liu.  An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations.  Communications in Computational Physics, 12(4):955–980, 2012.
  * [25] J. Hu, S. Jin, and Q. Li.  Asymptotic-preserving schemes for multiscale hyperbolic and kinetic equations.  Handbook of Numerical Analysis, 18:103–129, 2017.
  * [26] G.-S. Jiang and C.-W. Shu.  Efficient implementation of weighted ENO schemes.  Journal of Computational Physics, 126(1):202–228, 1996.
  * [27] S. Jin.  Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations.  SIAM Journal on Scientific Computing, 21(2):441–454, 1999.
  * [28] S. Klainerman and A. Majda.  Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids.  Communications on pure and applied Mathematics, 34(4):481–524, 1981.
  * [29] S. Klainerman and A. Majda.  Compressible and incompressible fluids.  Communications on Pure and Applied Mathematics, 35(5):629–651, 1982\.
  * [30] R. Klein.  Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics i: One-dimensional flow.  Journal of Computational Physics, 121(2):213–237, 1995.
  * [31] R. Klein.  An applied mathematical view of meteorological modelling.  In Applied Mathematics Entering the 21st century; Invited talks from the ICIAM 2003 Congress, volume 116, pages 227–269, 2004.
  * [32] R. Klein, E. Mikusky, and A. Owinoh.  Multiple scales asymptotics for atmospheric flows.  In Proceedings of the European Congress of Mathematics, pages 149–164. Springer, New York, 2011.
  * [33] A. Kurganov.  Finite-volume schemes for shallow-water equations.  Acta Numerica, 27:289–351, 2018.
  * [34] A. Kurganov and G. Petrova.  A second-order well-balanced positivity preserving central-upwind scheme for the saint-venant system.  Communications in Mathematical Sciences, 5(1):133–160, 2007.
  * [35] R. J. LeVeque.  Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm.  Journal of Computational Physics, 146(1):346–365, 1998.
  * [36] P. Li, W. S. Don, and Z. Gao.  High order well-balanced finite difference WENO interpolation-based schemes for shallow water equations.  Computers & Fluids, 201:104476, 2020.
  * [37] X. Liu.  A well-balanced asymptotic preserving scheme for the two-dimensional shallow water equations over irregular bottom topography.  SIAM Journal on Scientific Computing, 42(5):B1136–B1172, 2020.
  * [38] X. Liu.  A new well-balanced finite-volume scheme on unstructured triangular grids for two-dimensional two-layer shallow water flows with wet-dry fronts.  Journal of Computational Physics, 438:110380, 2021.
  * [39] X. Liu, A. Chertock, and A. Kurganov.  An asymptotic preserving scheme for the two-dimensional shallow water equations with coriolis forces.  Journal of Computational Physics, 391:259–279, 2019.
  * [40] A. Majda.  Introduction to PDEs and Waves for the Atmosphere and Ocean.  American Mathematical Society, 2003.
  * [41] F. Miczek, F. Röpke, and P. Edelmann.  A new numerical solver for flows at various Mach numbers.  Astronomy & Astrophysics, Vol. 576:A50, 2015.
  * [42] S. Noelle, G. Bispen, K. R. Arun, M. Lukáčová-Medvid́ová, and C.-D. Munz.  A weakly asymptotic preserving low mach number scheme for the euler equations of gas dynamics.  SIAM Journal on Scientific Computing, 36(6):B989–B1024, 2014.
  * [43] S. Noelle, N. Pankratz, G. Puppo, and J. R. Natvig.  Well-balanced finite volume schemes of arbitrary order of accuracy for shallow water flows.  Journal of Computational Physics, 213(2):474–499, 2006.
  * [44] S. Noelle, Y. Xing, and C.-W. Shu.  High-order well-balanced finite volume WENO schemes for shallow water equation with moving water.  Journal of Computational Physics, 226(1):29–58, 2007.
  * [45] L. Pareschi and G. Russo.  Implicit-Explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation.  Journal of Scientific Computing, 25(112):129–155, 2005.
  * [46] M. Ricchiuto.  An explicit residual based approach for shallow water flows.  Journal of Computational Physics, 280:306–344, 2015.
  * [47] M. Ricchiuto and A. Bollermann.  Stabilized residual distribution for shallow water simulations.  Journal of Computational Physics, 228(4):1071–1115, 2009.
  * [48] C.-W. Shu.  Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws.  In Advanced numerical approximation of nonlinear hyperbolic equations, pages 325–432. Springer, 1998.
  * [49] C.-W. Shu.  High order weighted essentially nonoscillatory schemes for convection dominated problems.  SIAM Review, 51(1):82–126, 2009.
  * [50] M. Tang.  Second order method for isentropic Euler equation in the low Mach number regime.  Kinetic and Related Models, 5(1):155–184, 2012.
  * [51] M. Tavelli and M. Dumbser.  A pressure-based semi-implicit space–time discontinuous Galerkin method on staggered unstructured meshes for the solution of the compressible Navier–Stokes equations at all Mach numbers.  Journal of Computational Physics, 341:341–376, 2017.
  * [52] A. Thomann, G. Puppo, and C. Klingenberg.  An all speed second order well-balanced IMEX relaxation scheme for the Euler equations with gravity.  Journal of Computational Physics, 420:109723, 2020.
  * [53] G. Tumolo, L. Bonaventura, and M. Restelli.  A semi-implicit, semi-Lagrangian, p-adaptive discontinuous Galerkin method for the shallow water equations.  Journal of Computational Physics, 232(1):46–67, 2013.
  * [54] E. Turkel.  Preconditioned methods for solving the incompressible and low speed compressible equations.  Journal of Computational Physics, 72:277–298, 1987.
  * [55] G. K. Vallis.  Atmospheric and Oceanic Fluid Dynamics.  Cambridge University Press, 2017.
  * [56] S. Vater and R. Klein.  A semi-implicit multiscale scheme for shallow water flows at low Froude number.  Communications in Applied Mathematics and Computational Science, 13(2):303–336, 2018.
  * [57] C. Viozat.  Implicit upwind schemes for low Mach number compressible flows.  PhD thesis, Inria, 1997.
  * [58] S. Vukovic and L. Sopta.  ENO and WENO schemes with the exact conservation property for one-dimensional shallow water equations.  Journal of Computational Physics, 179(2):593–621, 2002.
  * [59] Z. Wang, J. Zhu, and N. Zhao.  A new fifth-order finite difference well-balanced multi-resolution WENO scheme for solving shallow water equations.  Computers & Mathematics with Applications, 80(5):1387–1404, 2020\.
  * [60] X. Wen, W. S. Don, Z. Gao, and Y. Xing.  Entropy stable and well-balanced discontinuous Galerkin methods for the nonlinear shallow water equations.  Journal of Scientific Computing, 83(3):1–32, 2020.
  * [61] Y. Xing.  Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium.  Journal of Computational Physics, 257:536–553, 2014.
  * [62] Y. Xing and C.-W. Shu.  High order finite difference WENO schemes with the exact conservation property for the shallow water equations.  Journal of Computational Physics, 208(1):206–227, 2005.
  * [63] Y. Xing and C.-W. Shu.  A survey of high order schemes for the shallow water equations.  Journal of Mathematical Study, 47:221–249, 2014.
  * [64] Y. Xing, C.-W. Shu, and S. Noelle.  On the advantage of well-balanced schemes for moving-water equilibria of the shallow water equations.  Journal of Scientific Computing, 48(1):339–349, 2011.
  * [65] Y. Xing, X. Zhang, and C.-W. Shu.  Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations.  Advances in Water Resources, 33:1476–1493, 2010.
  * [66] J. Zeifang, J. Schütz, K. Kaiser, A. Beck, M. Lukáčová-Medvidová, and S. Noelle.  A novel full-Euler low Mach number IMEX splitting.  Communication in Computational Physics, 27:292–320, 2020.
  * [67] M. Zhang, W. Huang, and J. Qiu.  A high-order well-balanced positivity-preserving moving mesh DG method for the shallow water equations with non-flat bottom topography.  Journal of Scientific Computing, 87(3):1–43, 2021.