# High order asymptotic preserving and well-balanced schemes for the shallow
water equations with source terms

###### Abstract.

In this study, we investigate the Shallow Water Equations incorporating source
terms accounting for Manning friction and a non-flat bottom topology. Our
primary focus is on developing and validating numerical schemes that serve a
dual purpose: firstly, preserving all steady states within the model, and
secondly, maintaining the late-time asymptotic behavior of solutions, which is
governed by a diffusion equation and coincides with a long time and stiff
friction limit.

Our proposed approach draws inspiration from a penalization technique adopted
in [Boscarino et. al, SIAM Journal on Scientific Computing, 2014]. By
employing an additive implicit-explicit Runge-Kutta method, the scheme can
ensure a correct asymptotic behavior for the limiting diffusion equation,
without suffering from a parabolic-type time step restriction which often
afflicts multiscale problems in the diffusive limit. Numerical experiments are
performed to illustrate high order accuracy, asymptotic preserving, and
asymptotically accurate properties of the designed schemes.

###### Key words and phrases:

Shallow water equations; Manning friction; asymptotic preserving; well-
balanced; implicit-explicit Runge-Kutta; high order

Guanlan Huang

School of Mathematical Sciences, Xiamen University

Xiamen, Fujian 361005, PR China

<EMAIL_ADDRESS>
Sebastiano Boscarino

Department of Mathematics and Computer Science, University of Catania

Catania 95125, Italy

<EMAIL_ADDRESS>
Tao Xiong

School of Mathematical Sciences, Fujian Provincial Key Laboratory of
Mathematical Modeling

and High-Performance Scientific Computing, Xiamen University

Xiamen, Fujian 361005, PR China

<EMAIL_ADDRESS>

##  1\. Introduction

The Shallow Water Equations (SWEs) are a set of equations derived from the
more comprehensive Navier-Stokes equations, which govern the behavior of
fluids. This hyperbolic system finds widespread application in modeling fluid
dynamics in various natural environments, including channels, rivers, and
oceans. The study of these equations carries significant practical
implications, with numerous applications extending to crucial areas such as
the prediction of storm surge levels, tides, and alterations in coastlines
caused by events like hurricanes, among many others [46, 34, 49, 47, 28].

In recent years, the SWEs have become more and more important in atmospheric
flows, especially with different space and time scales. In this paper, we
consider the SWEs with Manning friction and a non-flat bottom topology, which
can be written as

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

where $h$ is the depth of the water layer, ${\bf u}={\bf q}/h$ is the flow
velocity, which are defined on a time-space domain $(t,{\bf
x})\in\mathbb{R}^{+}\times\Omega$; $g$ is the gravitational constant, $b({\bf
x})$ is the bottom topography which is independent of time; $\otimes$ denotes
the Kronecker product, and $\gamma{\bf q}$ models the bottom friction. Here
$\gamma=\dfrac{gk^{2}|{\bf q}|}{h^{\eta}}$ with $|{\bf q}|$ being the $L^{2}$
norm of ${\bf q}$, and $k$ is the Manning coefficient to determine the
intensity of the friction exerted by the bottom of the water. The larger the
$k$ is, the greater the friction exerted by the water. $\eta$ is a parameter
which is usually taken to be $7/3$.

To investigate the late-time behavior of solutions to the system (1.1), and
following the lines of [13], a behavior of $h$ and ${\bf q}$ in a long time
simulation and dominant friction is considered. Such an asymptotic behavior is
governed by a diffusive regime. For the convenience of analysis, we introduce
a small parameter $\varepsilon$ in order to scale the time $t$ and the
friction coefficient $k$, as follows:

(1.2)
$t\leftarrow\frac{\hat{t}}{\varepsilon}\,,\qquad{k\leftarrow\frac{\hat{k}}{\varepsilon}}.\,$

Then the system (1.1) can be rewritten as

(1.3)
$\left\\{\begin{array}[]{ll}\varepsilon\partial_{\hat{t}}h+\nabla\cdot{\bf
q}=0,\\\\[8.53581pt] \displaystyle\varepsilon\partial_{\hat{t}}{\bf
q}+\nabla\cdot\left(\dfrac{{\bf q}\otimes{\bf
q}}{h}\right)+\nabla\left(\dfrac{g}{2}h^{2}\right)=-gh\nabla
b-\frac{1}{\varepsilon^{2}}\widehat{\gamma}{\bf q},\end{array}\right.$

with $\widehat{\gamma}=\dfrac{g\hat{k}^{2}|{\bf q}|}{h^{\eta}}$. In [6] (and
the references therein), it was shown that, in the modeling of shallow water
problems when the friction is particularly strong and very much dominates over
the advection effects, the asymptotic regime, i.e. $\varepsilon\to 0$,
associated with (1.3) is governed by the nonlinear diffusion equation

(1.4) $h_{t}=\nabla\cdot\left(\frac{\sqrt{h}}{\kappa(h)}\frac{\nabla
h}{\sqrt{|\nabla h|}}\right).$

In the next section, we show that by the Chapman-Enskog expansions of the
variables ${\bf q}$ and $h$ in term of the parameter $\varepsilon$, we get the
limiting nonlinear diffusion equation (1.4) when $\varepsilon\to 0$.

Now, denoting

(1.5) $\widehat{{\bf q}}=\frac{{\bf q}}{\varepsilon},$

the system (1.3) becomes

(1.6)
$\left\\{\begin{array}[]{ll}\varepsilon\partial_{\hat{t}}h+\varepsilon\nabla\cdot\widehat{{\bf
q}}=0,\\\\[8.53581pt] \varepsilon^{2}\partial_{\hat{t}}\widehat{{\bf
q}}+\varepsilon^{2}\nabla\cdot\left(\dfrac{\widehat{{\bf
q}}\otimes\widehat{{\bf
q}}}{h}\right)+\nabla\left(\dfrac{g}{2}h^{2}\right)=-gh\nabla
b-\hat{\gamma}\widehat{{\bf q}}.\end{array}\right.$

Dividing the first sub-equation of system (1.6) by $\varepsilon$, the second
one by $\varepsilon^{2}$, and dropping the hats for clarity, then we can
obtain the non-dimensional system

(1.7) $\left\\{\begin{array}[]{ll}\partial_{t}h+\nabla\cdot{\bf
q}=0,\\\\[8.53581pt] \displaystyle\partial_{t}{\bf
q}+\nabla\cdot\left(\dfrac{{\bf q}\otimes{\bf
q}}{h}\right)+\frac{1}{\varepsilon^{2}}\nabla\left(\dfrac{g}{2}h^{2}\right)=-\frac{1}{\varepsilon^{2}}gh\nabla
b-\frac{1}{\varepsilon^{2}}\gamma{\bf q}.\end{array}\right.$

Notice that the nonlinear SWEs (1.7) belong to the family of hyperbolic
conservation laws with stiff source terms [6], i.e.,

(1.8) $U_{t}+\nabla\cdot F(U)=\frac{1}{\varepsilon^{2}}S(U),$

where $U=(h,{\bf q})^{T}$, $\nabla\cdot F(U)=\left(\nabla\cdot{\bf
q},\displaystyle\nabla\cdot\left(\frac{{\bf q}\otimes{\bf
q}}{h}\right)\right)^{T}$, and
$S(U)=-\left(0,\nabla\left(\dfrac{g}{2}h^{2}\right)+gh\nabla b+\gamma{\bf
q}\right)^{T}$.

Solving (1.8) numerically is a challenging task for many reasons. Firstly,
such system often admits non-trivial steady state solutions in which the
source term balances the effect of the flux gradients. There are many
interesting physical phenomena, which are small perturbations near these
steady states. A standard numerical scheme cannot capture this equilibrium
state or small perturbations near such equilibrium state, unless a much
refined mesh is adopted [46, 42, 45]. A significant scheme for the SWEs with
source terms is proposed by Bermudez and Vazquez [5], which balances the
gradient of pressure by source terms at the discrete level, so that the
equilibrium state can be maintained in coarse mesh subdivisions. Then many
schemes are designed to preserve the steady state [35, 53, 46, 49, 47, 48, 15,
22, 37, 39, 50, 45, 33], including finite difference, finite volume, and
discontinuous Galerkin methods. All the above works consider that the
solutions are not affected by the friction, while there are many situations
that the friction can not be ignored. The additional friction term makes the
well-balanced work more difficult, especially for situation that the friction
is dominant. When the system (1.7) reaches steady state, it becomes

(1.9) $\left\\{\begin{array}[]{ll}\nabla\cdot{\bf q}=0,\\\\[8.53581pt]
\displaystyle\nabla\cdot\left(\dfrac{{\bf q}\otimes{\bf
q}}{h}\right)+\frac{1}{\varepsilon^{2}}\nabla\left(\dfrac{g}{2}h^{2}\right)=-\frac{1}{\varepsilon^{2}}gh\nabla
b-\frac{1}{\varepsilon^{2}}\gamma{\bf q}.\end{array}\right.$

In case of a still water ${\bf q}=\mathbf{0}$, the system has a steady state
named “lake at rest” [52], that is

(1.10) $h({\bf x},t)+b({\bf x})=C,\qquad{\bf q}({\bf x},t)=\mathbf{0}.$

For a moving water steady state ${\bf q}\neq\mathbf{0}$, it makes the design
of a well-balanced scheme much more challenging. Very limited studies include
traditional time discretizations [36, 24], or semi-implicit time
discretizations [41, 44, 43, 17]. However, in most of these works, the stiff
friction is caused by the vanishing of water depth $h$. Another stiff friction
is caused by a large coefficient $k$ [13], which is considered in this work.
Besides, for simplicity, we only consider the still water steady state (1.10).

Secondly, a system as (1.8) with a small parameter $\varepsilon$ may relax to
a reduced system as $\varepsilon$ vanishes. Under suitable conditions [16],
the solutions of the system (1.8) will converge to the solutions of the
limiting system, see for example (1.4). It is not at all obvious that a
numerical scheme adopted for the solution of the system (1.8) will become (for
given space and time discretization parameters) a consistent discretization
for the limiting equation, as the small parameter $\varepsilon$ vanishes. When
this occurs, the scheme is said to be asymptotic preserving (AP). The AP
scheme was proposed by S. Jin [30] and we refer the readers to the review
papers [27, 31] for more work about AP schemes. Recently, Bulteau et. al. [13]
designed an AP scheme for system (1.7) by using an explicit time
discretization, combined with finite volume and a well-balanced modification
for space discretization. The scheme has been proven to be both AP and well-
balanced. However, due to an explicit time discretization, the scheme might be
subject to severe time step restrictions in order to maintain stability, as
the eigenvalues of the system being $\lambda_{1}={\bf
u}\cdot{\bf{n}}+c/\varepsilon$ and $\lambda_{2}={\bf
u}\cdot{\bf{n}}-c/\varepsilon$ along the normal direction ${\bf n}$, where
$c=\sqrt{gh}$ representing the speed of sound. Such eigenvalues lead to a time
step for an explicit scheme as

(1.11) $\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.

In the present paper, we will develop a numerical scheme to approximate the
solution of (1.7), which accurately captures the steady state solutions and is
able to restore the asymptotic regime, i.e., the scheme becomes a consistent
discretization for the limiting equation (1.4) as $\varepsilon\to 0$ (AP). Our
main focus is to design a class of AP implicit-explicit (IMEX) high order
finite difference weighted essentially non-oscillatory (WENO) schemes.
However, the AP property only guarantees the consistency of the scheme, and in
general does not guarantee the high order accuracy in the limit for
$\varepsilon\ll 1$, i.e. the order of accuracy may degrade. We will show that
under some assumptions, in the limit for $\varepsilon\to 0$, the scheme
maintains its order of temporal accuracy, and we say that the scheme is
asymptotically accurate (AA) for the limiting system when $\varepsilon\to 0$.
To achieve high order temporal accuracy, IMEX Runge-Kutta (RK) schemes [8, 7]
represent a powerful tool for the time discretization of stiff systems of the
form (1.8). However, directly applying IMEX RK schemes to the SWEs (1.7) are
not very efficient. In the diffusive regime, i.e. $\varepsilon\to 0$, such
schemes become consistent but explicit schemes for the diffusive limiting
equation (1.4). Therefore, these schemes suffer from a usual parabolic time
step restriction $\Delta t=O(\Delta x^{2})$, for given space and time
discretization parameters $\Delta x$ and $\Delta t$. In order to avoid this
parabolic time-step restriction, we introduce a new reformulation of the
problem (1.7) based on the addition of two opposite diffusive terms in the
first equation in (1.7), which has been introduced in [12, 9, 8, 51], and the
resulting scheme is consistent and relaxed to an implicit scheme in the
diffusive limit, therefore avoiding the typical parabolic restriction. This
strategy is very effective. Furthermore, we will show that, under this
reformulation of the system (1.7), the scheme converges to a high order
implicit scheme in the diffusive limit (1.4), namely the scheme is AA.

The rest of this paper is organized as follows. In Section 2, we will formally
derive the limiting equation for the SWEs with a non-flat topography and
Manning friction. An IMEX time discretization and finite difference WENO
reconstruction with well-balanced property will also be introduced. Formal
analyses of AP and AA properties of the scheme are provided in Section 3.
Numerical experiments are performed in Section 4, followed by concluding
remarks in Section 5.

##  2\. Diffusive limit and numerical scheme

In this section, we will first briefly review the diffusive limit for the SWEs
with strong friction and a bottom topography, which leads to a nonlinear
degenerate $p$-Laplacian equation. Then we will present an AP scheme, which
can well capture this nonlinear diffusive limit. This approach consists of
adding and subtracting the limiting nonlinear diffusive term in the water
depth equation, with different time discretizations and we will show that in
the diffusive limit the numerical scheme becomes an implicit time
discretization for the limiting system, which can avoid a parabolic time step
restriction. This technique has been developed in [12, 9, 8] and used in [38,
51].

For this system, since all terms with explicit and implicit time
discretizations are seperated, we employ an additive implicit-explicit (IMEX)
scheme, with a multi-stage RK method to achieve high order accuracy. In space,
we adopt a high order finite difference WENO reconstruction for convective
terms [40, 29], combined with central difference discretizations for second
order derivative terms. With a careful treatment in space, the well-balanced
property for a hydrostatic equilibrium state is also achieved. For our
proposed scheme, we can formally prove it is AP and AA.

###  2.1. Diffusive limit

Firstly, let us briefly review the limiting system for the SWEs with friction
and a non-flat bottom (1.7), under the diffusive scaling with (1.2) and (1.5).
We denote $H=h+b$ as the surface level of water, and assume the following
Chapman-Enskog expansions [13]

(2.1) $\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
q}({\bf x},t)={\bf q}_{0}({\bf x},t)+\varepsilon{\bf q}_{1}({\bf
x},t)+\varepsilon^{2}{\bf q}_{2}({\bf x},t)+\cdots,\end{array}\right.$

with

(2.2) $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.$

Here the bottom topography $b({\bf x})$ is independent of time. Substituting
(2.1) into (1.7), for the leading order terms, we have

(2.3a) $\displaystyle\partial_{t}h_{0}+\nabla\cdot{\bf q}_{0}=0,$ (2.3b)
$\displaystyle\nabla\left(\dfrac{g}{2}h_{0}^{2}\right)+gh_{0}\nabla
b=-\gamma_{0}{\bf q}_{0},$

with the nonlinear friction coefficient $\gamma_{0}=gk^{2}|{\bf
q}_{0}|/h_{0}^{\eta}$. Denoting $H_{0}({\bf x},t)=h_{0}({\bf x},t)+b({\bf
x})$, since

(2.4) $\nabla\left(\dfrac{g}{2}h_{0}^{2}\right)+gh_{0}\nabla b=gh_{0}\nabla
h_{0}+gh_{0}\nabla b=gh_{0}\nabla H_{0},$

with $\gamma_{0}=gk^{2}|{\bf q}_{0}|/h_{0}^{\eta}$, substituting (2.4) into
(2.3b) we obtain

(2.5) $gk^{2}|{\bf q}_{0}|^{2}/h_{0}^{\eta}=|\gamma_{0}{\bf
q}_{0}|=gh_{0}|\nabla H_{0}|,$

which yields $|{\bf q}_{0}|=\sqrt{h_{0}^{\eta+1}|\nabla H_{0}|/k^{2}}$, and
(2.3b) becomes

(2.6) ${\bf q}_{0}=-\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}.$

Substituting (2.6) into (2.3a), it leads to the following limiting equation

(2.7)
$\partial_{t}h_{0}=\nabla\cdot\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}\right).$

This limiting equation (2.7) is a non-stationary $p$-Laplacian equation with
$p=3/2$. $p$-Laplacian equations are degenerate parabolic equations, which
arise in many physical problems [20], and they have been studied both
theoretically and numerically in several works, such as [1, 26, 32, 19, 21,
18, 23, 4, 3], and references therein. Due to the limiting equation itself is
already nonlinear and degenerate, it makes the designing of an efficient
numerical scheme for (1.7) which can capture the diffusive limit (2.7) very
challenging. In the following, we will propose an efficient AP scheme, which
is high order both in space and in time and can well capture the corresponding
limit.

###  2.2. IMEX time discretization

We consider an IMEX time discretization for (1.7), which only uses implicit
time discretizations for stiff and diffusion terms, while convective terms are
discretized explicitly, so that we can achieve a large time step condition,
while avoiding to solve a very complicate nonlinear system. We start with
writing (1.7) in the following form of ordinary differential equations

(2.8) $U_{t}=F(U)+G(U),$

with $U=(h,{\bf q})^{T}$ and

(2.9) $F(U)=\begin{pmatrix}-\nabla\cdot{\bf
q}-\mu(\varepsilon)\nabla\cdot\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{\nabla
H}{\sqrt{|\nabla H|}}\right)\\\ \,\\\ -\nabla\cdot\left(\dfrac{{\bf
q}\otimes{\bf
q}}{h}\right)\end{pmatrix},\,G(U)=\begin{pmatrix}\mu(\varepsilon)\nabla\cdot\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{\nabla
H}{\sqrt{|\nabla H|}}\right)\\\ \,\\\
-\dfrac{1}{\varepsilon^{2}}\nabla\left(\dfrac{g}{2}h^{2}\right)-\dfrac{1}{\varepsilon^{2}}gh\nabla
b-\dfrac{1}{\varepsilon^{2}}\gamma{\bf q}\end{pmatrix},$

where we have added and substracted a term for the equation of water depth
$h$, based on the limiting equation (2.7). Here the weighted function
$\mu(\varepsilon)$ is such that $\mu:\mathbb{R}^{+}\to[0,1]$ and $\mu(0)=1$.
When $\varepsilon$ is not small, there is no need to add and subtract the
limiting term, therefore $\mu(\varepsilon)$ will be small in such a hyperbolic
regime, i.e., $\mu(\varepsilon)\approx 0$. The precise choice of
$\mu(\varepsilon)$ depends not only on $\varepsilon$ but also the mesh size or
the time step as discussed in [8], which will be specified later. Now let us
focus on the time discretization while keeping the space continuous. The
spatial discretizations will be discussed in the next section.

####  2.2.1. First order IMEX scheme

We begin with a first order IMEX scheme for (2.8) with (2.9). The scheme reads

(2.10a) $\displaystyle\dfrac{h^{n+1}-h^{n}}{\Delta t}+\nabla\cdot{\bf
q}^{n}+\mu(\varepsilon)\nabla\cdot\left(\sqrt{\dfrac{h^{\eta+1}}{k^{2}}}\dfrac{\nabla
H}{\sqrt{|\nabla
H|}}\right)^{n}=\mu(\varepsilon)\nabla\cdot\left(\sqrt{\dfrac{h^{\eta+1}}{k^{2}}}\frac{\nabla
H}{\sqrt{|\nabla H|}}\right)^{n+1},$ (2.10b) $\displaystyle\dfrac{{\bf
q}^{n+1}-{\bf q}^{n}}{\Delta t}+\nabla\cdot\left(\dfrac{{\bf q}\otimes{\bf
q}}{h}\right)^{n}+\dfrac{1}{\varepsilon^{2}}\nabla\left(\dfrac{g}{2}h^{2}\right)^{n+1}=-\dfrac{1}{\varepsilon^{2}}gh^{n+1}\nabla
b-\dfrac{1}{\varepsilon^{2}}\gamma^{n+1}{\bf q}^{n+1}.$

Here we take the weighted function to be
$\mu(\varepsilon)=e^{-\varepsilon^{2}/\Delta x}$ with $\Delta x$ being the
mesh size [12, 9]. This weighted function depends on $\varepsilon$ and $\Delta
x$, which ensures that as $\varepsilon$ goes to zero, $\mu(\varepsilon)$ will
approach $1$, while it approaches $0$ when $\varepsilon$ is of
$\mathcal{O}(1)$ and $\Delta x$ is small enough. In this case, as
$\varepsilon\rightarrow 0$, the second and third terms on the left hand side
cancel with each other, and the scheme (2.10a) formally converges to

(2.11) $\frac{h^{n+1}-h^{n}}{\Delta
t}-\nabla\cdot\left\\{\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{\nabla
H}{\sqrt{|\nabla H|}}\right\\}^{n+1}=0,$

which is an implicit time discretization for the limting equation (2.7), so
that a parabolic type time discretization namely $\Delta t=\mathcal{O}(\Delta
x^{2})$ can be avoided.

However, solving (2.10) we face two important issues. Let us rewrite (2.10a)
as

(2.12) $h^{n+1}=h_{\star}+\Delta
t\mu(\varepsilon)\nabla\cdot\left\\{\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{\nabla
H}{\sqrt{|\nabla H|}}\right\\}^{n+1},$

with

$h_{\star}=h^{n}-\Delta t\left\\{\nabla\cdot{\bf
q}^{n}+\mu(\varepsilon)\nabla\cdot\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{\nabla
H}{\sqrt{|\nabla H|}}\right)^{n}\right\\}.$

Since $H^{n+1}=h^{n+1}+b$, (2.12) is a fully nonlinear system of $h^{n+1}$,
especially the diffusion coefficient depends on $\sqrt{|\nabla H|}$ in the
denominator. The first issue is that when the water is at or close to rest, we
have $\nabla H=0$ or $\nabla H$ close to $0$. In this case, this nonlinear
diffusion term cannot be linearized by using the approach in [51]. Instead we
have to take $\nabla H$ in the numerator and denominator of (2.12) to be at
the same time level to solve a nonlinear system, in order to avoid such a
singularity in the denominator. We adopt a simple Picard iteration to solve
this nonlinear system as follows:

  * •

Starting from the initial guess $h^{n+1,0}=h^{n}$, at each iteration point, we
first compute $|\nabla H|$ (based on a numerical approximation), if $|\nabla
H|<\delta$ where $\delta$ is a small value (e.g., we take $\delta=10^{-9}$ in
our numerical tests, unless otherwise specified), we stop the iteration and
take $h^{n+1}=h^{n+1,0}=h^{n}$;

  * •

Otherwise, we continue the iteration by taking

(2.13) $h^{n+1,\ell+1}=h_{\star}+\Delta
t\mu(\varepsilon)\nabla\cdot\left\\{\frac{\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\right)^{n+1,\ell}}{\sqrt{|\nabla
H|^{n+1,\ell}}}\nabla H^{n+1,\ell+1}\right\\},\quad\ell\geq 0.$

The iteration is stopped when $\|h^{n+1,\ell+1}-h^{n+1,\ell}\|\leq\delta$.
Here $\delta$ is taken to be the same as above, and an $L^{1}$ norm is taken
in practice.

By taking this simple Picard iteration, as we can see, when initially the
solution is at a hydrodynamic equilibrium, that is $H=\text{Const.}$, due to
the iteration is stopped at the first step, such an equilibrium can be
preserved.

After obtaining $h^{n+1}$, the second issue is how to update ${\bf q}^{n+1}$
by (2.10b), since the friction coefficient $\gamma^{n+1}=\dfrac{gk^{2}|{\bf
q}^{n+1}|}{(h^{n+1})^{\eta}}$ also depends on ${\bf q}^{n+1}$. We first
rewrite (2.10b) as

(2.14) $\left(1+\frac{\Delta t}{\varepsilon^{2}}\,\frac{gk^{2}|{\bf
q}^{n+1}|}{(h^{n+1})^{\eta}}\right){\bf q}^{n+1}={\bf q}_{\star},$

here

(2.15) ${\bf q}_{\star}={\bf q}^{n}-\Delta t\nabla\cdot\left(\frac{{\bf
q}\otimes{\bf q}}{h}\right)^{n}-\frac{\Delta
t}{\varepsilon^{2}}\nabla\left(\frac{g}{2}h^{2}\right)^{n+1}-\frac{\Delta
t}{\varepsilon^{2}}gh^{n+1}\nabla b.$

Taking an $L^{2}$ norm on both sides of (2.14), it becomes a quadratic
equation for $|{\bf q}^{n+1}|$

(2.16) $\frac{\Delta t}{\varepsilon^{2}}\,\frac{gk^{2}}{(h^{n+1})^{\eta}}|{\bf
q}^{n+1}|^{2}+|{\bf q}^{n+1}|-|{\bf q}_{\star}|=0.$

(2.16) has two real roots, and the positive value is the reasonable one

$\displaystyle|{\bf
q}^{n+1}|=\frac{-\varepsilon^{2}+\sqrt{\varepsilon^{4}+\displaystyle\frac{4\Delta
tgk^{2}}{(h^{n+1})^{\eta}}|\varepsilon^{2}{\bf
q}_{\star}|}}{\displaystyle\frac{2\Delta tgk^{2}}{(h^{n+1})^{\eta}}}.$

Substituting it into (2.14), we update ${\bf q}^{n+1}$ by

(2.17) ${\bf q}^{n+1}=\frac{2\varepsilon^{2}{\bf
q}_{\star}}{\varepsilon^{2}+\sqrt{\varepsilon^{4}+\displaystyle\frac{4\Delta
tgk^{2}}{(h^{n+1})^{\eta}}|\varepsilon^{2}{\bf q}_{\star}|}}.$

However, we note that as $\varepsilon\rightarrow 0$, numerically in (2.17) we
face an issue of $\frac{0}{0}$, and in (2.15) a dividing by $\varepsilon^{2}$.
Both will lead to numerical instabilities. To avoid these two stiff issues,
from (2.15), we compute $\varepsilon^{2}{\bf q}_{\star}$ into one by

$\varepsilon^{2}{\bf q}_{\star}=\varepsilon^{2}{\bf
q}^{n}-\varepsilon^{2}\Delta t\nabla\cdot\left(\frac{{\bf q}\otimes{\bf
q}}{h}\right)^{n}-\Delta t\nabla\left(\frac{g}{2}h^{2}\right)^{n+1}-\Delta
tgh^{n+1}\nabla b.$

If the solution is at a hydrostatic equilibrium, namely
$\nabla\left(\dfrac{1}{2}h^{2}\right)=-h\nabla b$ and ${\bf q}=\mathbf{0}$, we
would expect ${\bf q}$ always to be $\mathbf{0}$. However, when $\varepsilon$
is small, due to round-off errors, the updating of ${\bf q}^{n+1}$ by (2.17)
is unstable. In this case, if $|\varepsilon^{2}{\bf q}_{\star}|<\eta$ and
$|{\bf q}^{n}|<\eta$ where $\eta$ is taken to be $10^{-12}$ at the round-off
error level, we set ${\bf q}^{n+1}=\mathbf{0}$.

####  2.2.2. High order IMEX scheme

The first order IMEX scheme (2.10) can be generalized to high order with an
additive IMEX RK method, such procedure has been introduced in [8]. For the
additive system (2.8), an additive IMEX RK time discretization is
characterized by the double Butcher $tableau$ [14],

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

where $\tilde{A}=(\tilde{a}_{ij})$ is an $s\times s$ lower triangular matrix,
i.e. $\tilde{a}_{ij}=0$ for $j\geq i$, leading to an explicit scheme, while
$A=({a}_{ij})$ is an $s\times s$ matrix for an implicit scheme. A diagonally
implicit RK (DIRK) scheme with $a_{ij}=0$ for $j>i$ is used, which is simple
and efficient for implicit discretization. Other vectors are
$\tilde{\mathbf{c}}=(\tilde{c}_{1},...,\tilde{c}_{s})^{T}$,
$\tilde{\mathbf{b}}=(\tilde{b}_{1},...,\tilde{b}_{s})^{T}$, and
$\mathbf{c}=(c_{1},...,c_{s})^{T}$, $\mathbf{b}=(b_{1},...,b_{s})^{T}$, where
$\tilde{\mathbf{c}}$ and $\mathbf{c}$ satisfy the relation
$\tilde{c}_{i}=\sum_{j=1}^{i-1}\tilde{a}_{ij}$, $c_{i}=\sum_{j=1}^{i}a_{ij}$
for $i=1,...,s$.

A sketched procedure of the additive IMEX RK scheme is given as follows.
Starting from $U^{n}$, for inner stages $i=1,\ldots,s$:

  * •

Compute the solution $U_{\star}^{(i)}$ from previous known stages

(2.19) $U_{\star}^{(i)}=U^{n}+\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}F(U^{(j)})+\Delta
t\sum_{j=1}^{i-1}a_{ij}G(U^{(j)}).$

  * •

Update $U^{(i)}$ with

(2.20) $U^{(i)}=U_{\star}^{(i)}+\Delta ta_{ii}G(U^{(i)}).$

After obtaining all inner stage values, we update the solution $U^{n+1}$ from

(2.21) $U^{n+1}=U^{n}+\Delta t\sum_{i=1}^{s}\tilde{b}_{i}F(U^{(i)})+\Delta
t\sum_{i=1}^{s}b_{i}G(U^{(i)}).$

Applying (2.19)-(2.21) to (2.8) with (2.9), in components, for
$U_{\star}^{(i)}$ in (2.19), they are

(2.22a) $\displaystyle h_{\star}^{(i)}=\,$ $\displaystyle h^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\left\\{\nabla\cdot{\bf
q}^{(j)}+\mu(\varepsilon)\nabla\cdot\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{\nabla
H}{\sqrt{|\nabla H|}}\right)^{(j)}\right\\}$ $\displaystyle+\Delta
t\sum_{j=1}^{i-1}a_{ij}\mu(\varepsilon)\nabla\cdot\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{\nabla
H}{\sqrt{|\nabla H|}}\right)^{(j)},$ (2.22b) $\displaystyle{\bf
q}^{(i)}_{\star}=\,$ $\displaystyle{\bf q}^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\nabla\cdot\left(\frac{{\bf q}\otimes{\bf
q}}{h}\right)^{(j)}-\Delta
t\sum_{j=1}^{i-1}a_{ij}\left\\{\frac{1}{\varepsilon^{2}}\nabla\left(\frac{g}{2}h^{2}\right)+\frac{1}{\varepsilon^{2}}gh\nabla
b+\frac{1}{\varepsilon^{2}}\gamma{\bf q}\right\\}^{(j)},$

and for $U^{(i)}$ in (2.20), they are

(2.23a) $\displaystyle h^{(i)}=h_{\star}^{(i)}+a_{ii}\Delta
t\mu(\varepsilon)\nabla\cdot\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{\nabla
H}{\sqrt{|\nabla H|}}\right)^{(i)},$ (2.23b) $\displaystyle{\bf q}^{(i)}={\bf
q}_{\star}^{(i)}-a_{ii}\frac{\Delta
t}{\varepsilon^{2}}\left\\{\nabla\left(\frac{g}{2}h^{2}\right)+gh\nabla
b+\gamma{\bf q}\right\\}^{(i)}.$

The nonlinear system (2.23a) is also solved for $h^{(i)}$ by using a simple
Picard iteration, which is almost the same as the first order scheme. The
updating ${\bf q}^{(i)}$ is similar to (2.17), here it reads

(2.24) ${\bf q}^{(i)}=\frac{2\varepsilon^{2}{\bf
q}_{\star\star}^{(i)}}{\varepsilon^{2}+\sqrt{\varepsilon^{4}+\frac{4a_{ii}\Delta
tgk^{2}}{(h^{(i)})^{\eta}}|\varepsilon^{2}{\bf q}_{\star\star}^{(i)}|}},$

with

$\varepsilon^{2}{\bf q}_{\star\star}^{(i)}=\varepsilon^{2}{\bf
q}_{\star}^{(i)}-a_{ii}\Delta
t\left\\{\nabla\left(\frac{g}{2}h^{2}\right)+gh\nabla b\right\\}^{(i)}.$

For the high order IMEX RK scheme, as discussed in [10, 11], a stability
property is needed to ensure that the scheme is both AP and AA. Specifically,
we require that the IMEX RK scheme is _globally stiffly accurate_ (GSA) with
$\tilde{\mathbf{b}}^{T}={\bf e}_{s}^{T}\tilde{A}$ for the explicit part of the
scheme and $\mathbf{b}^{T}={\bf e}_{s}^{T}A$ for the implicit one where ${\bf
e}_{s}=(0,\cdots,0,1)^{T}$ [25]. Moreover, when the IMEX RK scheme is GSA, the
numerical solution coincides with the last stage of the method, i.e.,
$U^{n+1}=U^{(s)}$. From now on, we only consider GSA IMEX RK schemes.

###  2.3. High order spatial discretizations

For spatial approximations of (2.9), we follow [28] by using a finite
difference WENO reconstruction for convective terms, while a high order
compact difference method for second order diffusion terms. To preserve a
hydrostatic equilibrium state, a well-balanced reconstruction [46] between the
gradient of pressure and bottom topography is adopted. Unlike the work [28],
we do not need to form an elliptic equation for the system (1.7) in the
diffusive scaling.

Without loss of generality, we take the 2D problem of (2.9) as an example to
describe our spatial approximations. 1D and higher dimensions can be
discretized similarly. Denoting $U=(h,hu,hv)^{T}$, the fluxes in (2.9) can be
written as

$\displaystyle F(U)=F_{x}(U)+F_{y}(U),$ $\displaystyle
G(U)=G_{x}(U)+G_{y}(U),$

with

$\displaystyle
F_{x}(U)=\left(-(hu)_{x}-\mu(\varepsilon)\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{H_{x}}{\sqrt{H_{x}+H_{y}}}\right)_{x},\,\,-(hu^{2})_{x},\,\,-(huv)_{x}\right)^{T},$
$\displaystyle
F_{y}(U)=\left(-(hv)_{y}-\mu(\varepsilon)\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{H_{y}}{\sqrt{H_{x}+H_{y}}}\right)_{y},\,\,-(huv)_{y},\,\,-(hv^{2})_{y}\right)^{T},$
$\displaystyle
G_{x}(U)=\left(\mu(\varepsilon)\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{H_{x}}{\sqrt{H_{x}+H_{y}}}\right)_{x},\,\,-\frac{1}{\varepsilon^{2}}\left(\frac{g}{2}h^{2}\right)_{x}-\frac{1}{\varepsilon^{2}}ghb_{x}-\frac{1}{\varepsilon^{2}}\frac{gk^{2}\sqrt{(hu)^{2}+(hv)^{2}}}{h^{\eta}}(hu),\,\,0\right)^{T},$
$\displaystyle
G_{y}(U)=\left(\mu(\varepsilon)\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{H_{y}}{\sqrt{H_{x}+H_{y}}}\right)_{y},\,\,0,\,\,-\frac{1}{\varepsilon^{2}}\left(\frac{g}{2}h^{2}\right)_{y}-\frac{1}{\varepsilon^{2}}ghb_{y}-\frac{1}{\varepsilon^{2}}\frac{gk^{2}\sqrt{(hu)^{2}+(hv)^{2}}}{h^{\eta}}(hu)\right)^{T}.$

For simplicity, we consider a rectangular domain $\Omega=[a,b]\times[c,d]$,
with a uniform Cartesian partition. Let $\Delta x=\frac{b-a}{N_{x}}$ and
$\Delta y=\frac{d-c}{N_{y}}$ be the mesh sizes along the $x$ and $y$
directions respectively, where $N_{x}$ and $N_{y}$ are the number of cells
correspondingly. The cell is denoted as
$I_{ij}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]$
with $i=1,2,\ldots,N_{x}$ and $j=1,2,\ldots,N_{y}$, where
$a=x_{\frac{1}{2}}<x_{\frac{3}{2}}<\cdots<x_{N_{x}+\frac{1}{2}}=b$ and
$c=y_{\frac{1}{2}}<y_{\frac{3}{2}}<\cdots<y_{N_{y}+\frac{1}{2}}=d$. Now we set
the grid point $(x_{i},y_{j})$ as
$x_{i}=\frac{1}{2}(x_{i-\frac{1}{2}}+x_{i+\frac{1}{2}})$ and
$y_{j}=\frac{1}{2}(y_{j-\frac{1}{2}}+y_{j+\frac{1}{2}})$, located at each cell
center. We denote $\omega_{i,j}$ as the value at the point $(x_{i},y_{j})$,
where $\omega$ is $h$ or ${\bf q}$. Correspondingly
$\omega^{\pm}_{i+\frac{1}{2},j}$, $\omega^{\pm}_{i,j+\frac{1}{2}}$ represent
cell interface values at $(x_{i+\frac{1}{2}},y_{j})$ and
$(x_{i},y_{j+\frac{1}{2}})$ from the left and right limits, respectively. A
high order finite difference WENO reconstruction will be used for those
convective terms, which is done in a dimension-by-dimension way for the 2D
problems:

  * •

WENO reconstruction with a Lax-Friedrichs flux splitting.
Here we briefly review a finite difference WENO reconstruction with a Lax-
Friedrichs flux splitting for the convective terms $\nabla\cdot{\bf q}$ and
$\nabla\cdot\left(\frac{{\bf q}\otimes{\bf q}}{h}\right)$, we refer to [29,
40] for more details. Taking $\nabla\cdot{\bf q}$ as an example, it is in the
form

(2.25) $\nabla\cdot{\bf q}=(hu)_{x}+(hv)_{y},$

and we approximate them with a $k$th order finite difference reconstruction,
namely

(2.26a) $\displaystyle(hu)_{x}\big{|}_{(x_{i},y_{j})}=\frac{1}{\Delta
x}\left\\{\widehat{(hu)}_{i+\frac{1}{2},j}-\widehat{(hu)}_{i-\frac{1}{2},j}\right\\}+\mathcal{O}(\Delta
x^{k}),$ (2.26b) $\displaystyle(hv)_{y}\big{|}_{(x_{i},y_{j})}=\frac{1}{\Delta
y}\left\\{\widehat{(hv)}_{i,j+\frac{1}{2}}-\widehat{(hv)}_{i,j-\frac{1}{2}}\right\\}+\mathcal{O}(\Delta
y^{k}).$

${\widehat{(hu)}_{i+\frac{1}{2},j}}$ and $\widehat{(hv)}_{i,j+\frac{1}{2}}$
are obtained by a high order finite difference WENO reconstruction. For
example, ${\widehat{(hu)}_{i+\frac{1}{2},j}}$ is decomposed into the following
two parts

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

$\widehat{(hu)}_{i+\frac{1}{2},j}^{+}$ and
$\widehat{(hu)}_{i+\frac{1}{2},j}^{-}$ are reconstructed by using upwind and
downwind information based on a Lax-Friedrichs flux splitting, with

(2.28) $(hu)_{i,j}^{\pm}=\frac{1}{2}\left\\{(hu)_{i,j}\pm\Lambda
H\right\\},\quad\Lambda=\max\left\\{|{\bf
u}|+\min\\{1,\frac{1}{\varepsilon^{2}}c\\}\right\\}.$

For example, for a fifth order WENO reconstruction, the upwind flux
$\widehat{(hu)}_{i+\frac{1}{2},j}^{+}$ is obtained from
$\\{(hu)_{i-2,j}^{+},(hu)_{i-1,j}^{+},(hu)_{i,j}^{+},(hu)_{i+1,j}^{+},(hu)_{i+2,j}^{+}\\}$,
while the downwind flux $\widehat{(hu)}_{i+\frac{1}{2},j}^{-}$ is obtained
based on
$\\{(hu)_{i-1,j}^{-},(hu)_{i,j}^{-},(hu)_{i+1,j}^{-},(hu)_{i+2,j}^{-},(hu)_{i+3,j}^{-}\\}$,
respectively. Note that here we take the surface level $H$ in (2.28) for
numerical viscosity, which in the hydrostatic state (1.10) will not introduce
any numerical dissipation, so that can help to preserve such an equilibrium
steady state. $(hv)_{i,j}^{\pm}$ can be obtained similarly, and it is also
similar for the flux difference approximation of $\nabla\cdot\left(\frac{{\bf
q}\otimes{\bf q}}{h}\right)$. We denote such a WENO reconstruction as
$\nabla_{LW}$.

  * •

Well-balanced WENO reconstrction.
If initially the water is at a hydrostatic state (1.10), such an equilibrium
will be preserved for later times. Numerically it is also very important to
keep such an equilibrium state. Many well-balanced methods have been proposed
[34, 53, 35, 46, 45]. Here we follow [46] to obtain a well-balanced WENO
reconstruction. For the system, in order to preserve ${\bf q}=\mathbf{0}$ when
$H$ is a constant, it is important to preserve

$\nabla\left(\frac{g}{2}h^{2}\right)=-gh\nabla b.$

The main idea is to split the right side as

$\nabla\left(\frac{g}{2}h^{2}\right)=-gh\nabla b=-gH\nabla
b+\nabla\left(\frac{g}{2}b^{2}\right),$

and apply the same WENO reconstruction (with exactly the same nonlinear
weights) of $\nabla\left(\frac{g}{2}h^{2}\right)$ to $\nabla b$ and
$\nabla\left(\frac{g}{2}b^{2}\right)$, so that

$\nabla\left(\frac{g}{2}h^{2}\right)+gH\nabla
b-\nabla\left(\frac{g}{2}b^{2}\right)=\nabla\left(\frac{g}{2}H^{2}\right)$

also holds at the discrete level, when $H$ is a constant. Such a WENO
reconstruction is denoted as $\nabla_{W}^{WB}$.

  * •

Compact central difference discretization for diffusive terms
For the diffusive terms in (2.9), a compact central difference discretization
will be applied [10, 11, 28]. For example, a fourth order central difference
for $(a(x,y)H_{x})_{x}$ at $(x_{i},y_{j})$ takes the following form

(2.29) $\displaystyle(a(x,y)H_{x})_{x}\Big{|}_{(x_{i},y_{j})}=$
$\displaystyle\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 H}^{T}_{i,j}$
$\displaystyle+\mathcal{O}(\Delta x^{4}),$

with the two vectors being

$\displaystyle{\bf a}_{i,j}$
$\displaystyle=(a_{i-2,j},a_{i-1,j},a_{i,j},a_{i+1,j},a_{i+2,j}),$
$\displaystyle{\bf H}_{i,j}$
$\displaystyle=(H_{i-2,j},H_{i-1,j},H_{i,j},H_{i+1,j},H_{i+2,j}),$

similarly for $(a(x,y)H_{y})_{y}$ along the $y$ direction. The function
$a(x,y)$ represents the coefficient
$\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{1}{\sqrt{|\nabla H|}}$. For the gradient
$|\nabla H|$ in the coefficient of the diffusive terms, as no upwind
information is available, so a WENO reconstruction without numerical viscosity
(namely $\Lambda=0$ as in (2.28)) is used, which we denote as $|\nabla_{W}H|$.
We note that as we observed numerically a WENO reconstruction of $|\nabla H|$
is important to control numerical oscillations, especially when the solution
is not smooth.

###  2.4. Algorithm flowchart

Now we will summarize our proposed scheme, with high order finite difference
WENO recontructions and compact central difference discretizations, coupled
with high order IMEX RK schemes, which is AP and well-balanced for the SWEs
with a nonflat bottom topology and Manning friction (1.7). The algorithm
flowchart is as follows.

  1. (1)

Start from the time level $t^{n}$, with $U^{n}$, we first compute
$U^{(i)}_{\star}$ from (2.19):

(2.30a) $\displaystyle h_{\star}^{(i)}=\,$ $\displaystyle h^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\left\\{\nabla_{LW}\cdot{\bf
q}^{(j)}+\mu(\varepsilon)\nabla\cdot\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{\nabla
H}{\sqrt{|\nabla_{W}H|}}\right)^{(j)}\right\\}$ $\displaystyle+\Delta
t\sum_{j=1}^{i-1}a_{ij}\mu(\varepsilon)\nabla\cdot\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{\nabla
H}{\sqrt{|\nabla_{W}H|}}\right)^{(j)},$ (2.30b) $\displaystyle{\bf
q}^{(i)}_{\star}=\,$ $\displaystyle{\bf q}^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\nabla_{LW}\cdot\left(\frac{{\bf q}\otimes{\bf
q}}{h}\right)^{(j)}-\Delta
t\sum_{j=1}^{i-1}a_{ij}\left\\{\frac{1}{\varepsilon^{2}}\left[\nabla_{W}^{WB}\left(\frac{g}{2}h^{2}\right)\right.\right.$
$\displaystyle\left.\left.+g(h+b)\nabla_{W}^{WB}b-\nabla_{W}^{WB}\left(\frac{g}{2}b^{2}\right)\right]+\frac{1}{\varepsilon^{2}}\frac{gk^{2}|{\bf
q}|}{h^{\eta}}{\bf q}\right\\}^{(j)}.$

  2. (2)

Then we solve $U_{I}^{(i)}$ from (2.20):

    1. (a)

$U_{I}^{(i)}$ can be expressed as

(2.31a) $\displaystyle h^{(i)}=\,$ $\displaystyle
h_{\star}^{(i)}+a_{ii}\mu(\varepsilon)\Delta
t\nabla\cdot\left(\sqrt{\frac{h^{\eta+1}}{k^{2}}}\frac{\nabla
H}{\sqrt{|\nabla_{W}H|}}\right)^{(i)},$ (2.31b) $\displaystyle{\bf
q}^{(i)}=\,$ $\displaystyle{\bf q}_{\star}^{(i)}-\Delta
ta_{ii}\left\\{\frac{1}{\varepsilon^{2}}\left[\nabla_{W}^{WB}\left(\frac{g}{2}h^{2}\right)+g(h+b)\nabla_{W}^{WB}b\right.\right.$
$\displaystyle\left.\left.-\nabla_{W}^{WB}\left(\frac{g}{2}b^{2}\right)\right]+\frac{1}{\varepsilon^{2}}\frac{gk^{2}|{\bf
q}|}{(h^{\eta})}{\bf q}\right\\}^{(i)}.$

    2. (b)

For the system (2.31), we first obtain $h^{(i)}$ from (2.31a) by using a
Picard iteration as the first order IMEX scheme, that is, we set
$h^{(i),0}=h^{(i-1)}$ and compute $|\nabla_{W}H^{(i),0}|$ with a WENO
reconstruction $\nabla_{W}$. If it is less than $\delta$, we set
$h^{(i)}=h^{(i),0}$ and stop the iteration, otherwise, we continue the
iteration by

(2.32) $h^{(i),\ell+1}=h_{\star}^{(i)}-a_{ii}\mu(\varepsilon)\Delta
t\nabla\cdot\left\\{\left(\frac{\sqrt{\frac{h^{\eta+1}}{k^{2}}}}{\sqrt{|\nabla_{W}H|}}\right)^{(i),\ell}\nabla
H^{(i),\ell+1}\right\\},\quad\ell\geq 0.$

The elliptic operator
$\nabla\cdot\left\\{\left(\frac{\sqrt{\frac{h^{\eta+1}}{k^{2}}}}{\sqrt{|\nabla_{W}H|}}\right)^{(i),\ell}\nabla
H^{(i),\ell+1}\right\\}$ in (2.32) is discreted by a high order central
difference discretization as described in (2.29). The iteration is stoped when
$||h^{(i),\ell+1}-h^{(i),\ell}||\leq\delta$.

    3. (c)

After obtaining $h^{(i)}$, we update ${\bf q}^{(i)}$ from

(2.33) ${\bf q}^{(i)}=\frac{2\varepsilon^{2}{\bf
q}_{\star\star}^{(i)}}{\varepsilon^{2}+\sqrt{\varepsilon^{4}+\displaystyle\frac{4a_{ii}\Delta
tgk^{2}}{(h^{(i)})^{\eta}}|\varepsilon^{2}{\bf q}_{\star\star}^{(i)}|}},$

with

(2.34) $\varepsilon^{2}{\bf q}_{\star\star}^{(i)}=\varepsilon^{2}{\bf
q}_{\star}^{(i)}-a_{ii}\Delta
t\left\\{\nabla_{W}^{WB}\left(\frac{g}{2}h^{2}\right)+g(h+b)\nabla_{W}^{WB}b-\nabla_{W}^{WB}\left(\frac{g}{2}b^{2}\right)\right\\}^{(i)}.$

Similarly, when $|\varepsilon^{2}{\bf q}_{\star\star}^{(i)}|$ and $|{\bf
q}^{(i-1)}|$ both are less than $\delta$, we set ${\bf q}^{(i)}=\mathbf{0}$.

  3. (3)

Finally, we update the numerical solution $U^{n+1}=U^{(s)}$ with the
assumption of IMEX RK schemes being GSA.

##  3\. Asymptotic preserving (AP) and Asymptotically Accurate (AA)

In this section, we will formally show that the first order IMEX scheme (2.10)
is AP, and the high order IMEX RK scheme (2.19)-(2.21) is AA. For both
properties, we only consider time discretizations, while keeping space
continuous. In accordance with (2.1), we adopt the following Chapman-Enskog
expansions:

(3.1) $h^{n}({\bf x})=h^{n}_{0}({\bf x})+\mathcal{O}(\varepsilon),\qquad{\bf
q}^{n}({\bf x})={\bf q}^{n}_{0}({\bf x})+\mathcal{O}(\varepsilon),$

with $H^{n}_{0}({\bf x})=h^{n}_{0}({\bf x})+b({\bf x})$, and

(3.2) ${\bf q}^{n}_{0}=-\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}\right)^{n}.$

Here $h^{n}({\bf x})=h({\bf x},t^{n}),\,{\bf q}^{n}({\bf x})={\bf q}({\bf
x},t^{n})$, and we say that the initial conditions are well-prepared if the
condition (3.1) and (3.2) are satisfied at $n=0$.

###  3.1. AP property.

For the first order IMEX scheme (2.10), we have

######  Theorem 3.1.

As $\varepsilon\rightarrow 0$ the first order IMEX scheme (2.10) is asymptotic
preserving as long as the initial conditions are well-prepared.

###### Proof.

Such a proof is straightforward by using the mathematical induction. As
$\varepsilon$ goes to zero, the system (2.10) becomes

(3.3a) $\displaystyle\frac{h_{0}^{n+1}-h_{0}^{n}}{\Delta t}+\nabla\cdot{\bf
q}^{n}_{0}+\mu(0)\nabla\cdot\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla
H_{0}|}}\right)^{n}=\mu(0)\nabla\cdot\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}\right)^{n+1},$ (3.3b)
$\displaystyle\nabla\left(\frac{g}{2}h_{0}^{2}\right)^{n+1}=-gh_{0}^{n+1}\nabla
b-\frac{gk^{2}|{\bf q}_{0}^{n+1}|}{h_{0}^{\eta}}{\bf q}_{0}^{n+1}.$

From (3.3b), following (2.4)-(2.6), we obtain

(3.4) ${\bf
q}_{0}^{n+1}=-\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}\right)^{n+1}.$

Substituting the relation of (3.4) at time level $t^{n}$ into (3.3a), since
$\mu(0)=1$, as $\varepsilon\to 0$ we get

(3.5) $\frac{h_{0}^{n+1}-h_{0}^{n}}{\Delta
t}=\nabla\cdot\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}\right)^{n+1},$

which is a consistent implicit time discretization to the limiting equation
(2.7). Starting from well-prepared initial conditions (3.1) at $n=0$, by the
mathematical induction, we can see the AP property holds for all later time
$t^{n}$. ∎

###  3.2. AA property.

Similarly as the AP property for the first order IMEX scheme (2.10), here we
focus on the AA analysis for high order IMEX RK schemes (2.19)-(2.21). In what
follows, we recognize that the GSA condition is crucial to guarantee the AA
property.

######  Theorem 3.2.

Consider applying an IMEX RK scheme (2.19)-(2.21) of order $p$. Assume that
the IMEX RK method is GSA and the the initial conditions are well-prepared.
Let us denote by ${\bf U}^{1}({\bf x};\varepsilon)=(h^{1}({\bf
x};\varepsilon),{\bf q}^{1}({\bf x};\varepsilon))^{T}$ the numerical solution
after one time step and ${\bf U}^{exa}({\bf x},t)=(h^{exa}({\bf x},t),{\bf
q}^{exa}({\bf x},t))^{T}$ the exact solution of the limiting equation (2.7).
Then after one step, the error between the numerical and exact solution
satisfies the following estimate

(3.6) $\lim_{\varepsilon\to 0}\,{\bf U}^{1}({\bf x};\varepsilon)={\bf
U}^{exc}({\bf x},\Delta t)+\mathcal{O}(\Delta t^{p+1}),$

namely, the scheme is AA.

###### Proof.

We prove it by using the mathematical induction for time stages. If we assume
(3.1) and (3.2) hold at $t^{n}$, from the AP property of the first order IMEX
scheme (2.10), we can get (3.4) and (3.5) hold for a time step $a_{11}\Delta
t$ at the first stage. Then if we assume this formal convergence holds for the
$(i-1)$-th stage, we will prove it also holds for the $i$-th stage. We have:

  1. (1)

for $U_{\star}^{(i)}$ in (2.22), as $\varepsilon\rightarrow 0$, it becomes

$\displaystyle h_{\star}^{(i)}=\,$ $\displaystyle h_{0}^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\left\\{\nabla\cdot{\bf
q}_{0}^{(j)}+\mu(0)\nabla\cdot\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}\right)^{(j)}\right\\},$ $\displaystyle+\Delta
t\sum_{j=1}^{i-1}a_{ij}\mu(0)\nabla\cdot\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}\right)^{(j)},$ $\displaystyle{\bf
q}^{(i)}_{\star}=\,$ $\displaystyle{\bf q}_{0}^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\nabla\cdot\left(\frac{{\bf q}_{0}\otimes{\bf
q}_{0}}{h_{0}}\right)^{(j)}.$

Since

${\bf q}_{0}^{(j)}=-\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}\right)^{(j)},\quad 0\leq j\leq i-1,$

and $\mu(0)=1$ as $\varepsilon\to 0$, $U^{(i)}_{\star}$ becomes

$\displaystyle h_{\star}^{(i)}=\,$ $\displaystyle h_{0}^{n}+\Delta
t\sum_{j=1}^{i-1}a_{ij}\nabla\cdot\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}\right)^{(j)},$ $\displaystyle{\bf
q}^{(i)}_{\star}=\,$ $\displaystyle{\bf q}_{0}^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\nabla\cdot\left(\frac{{\bf q}_{0}\otimes{\bf
q}_{0}}{h_{0}}\right)^{(j)}.$

  2. (2)

for $U^{(i)}$ in (2.23), as $\varepsilon\rightarrow 0$, first from (2.23b) we
get

$\nabla\left(\frac{g}{2}(h_{0}^{(i)})^{2}\right)+gh_{0}^{(i)}\nabla
b+\gamma^{(i)}_{0}{\bf q}_{0}^{(i)}=0,$

with $\gamma^{(i)}_{0}=\dfrac{gk^{2}|{\bf
q}_{0}^{(i)}|}{(h_{0}^{\eta})^{(i)}}$, so that

${\bf q}_{0}^{(i)}=-\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}\right)^{(i)},$

and (2.23a) becomes

$h_{0}^{(i)}=h_{\star}^{(i)}+\Delta
ta_{ii}\nabla\cdot\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}\right)^{(i)}=h_{0}^{n}+\Delta
t\sum_{j=1}^{i}a_{ij}\nabla\cdot\left(\sqrt{\frac{h_{0}^{\eta+1}}{k^{2}}}\frac{\nabla
H_{0}}{\sqrt{|\nabla H_{0}|}}\right)^{(j)}.$

Therefore, the formal convergence also holds for the stage $i$.

With the inner stage $i$ up to $s$, and given that the IMEX RK scheme is GSA,
we update the solution $U^{n+1}$ as $U^{n+1}=U^{(s)}$. Therefore, as
$\varepsilon\rightarrow 0$, the scheme becomes a $p$th order implicit scheme
for the limiting equation (2.7) using the matrix $A$ in (2.18). This means
that the error after one time step satisfies (3.6), confirming that the scheme
is AA. ∎

##  4\. Numerical tests

In this section, we will perform some numerical experiments for the SWEs with
Manning friction and a non-flat bottom topology, with our proposed scheme. A
third order GSA IMEX RK scheme from [2] is adopted, with the double Butcher
tableau given by

$\displaystyle\begin{array}[]{c|ccccc|ccccc}0&0&0&0&0&0&0&0&0&0&0\\\\[8.53581pt]
\frac{1}{2}&\frac{1}{2}&0&0&0&0&0&\frac{1}{2}&0&0&0\\\\[8.53581pt]
\frac{2}{3}&\frac{11}{18}&\frac{1}{18}&0&0&0&0&\frac{1}{6}&\frac{1}{2}&0&0\\\\[8.53581pt]
\frac{1}{2}&\frac{5}{6}&-\frac{5}{6}&\frac{1}{2}&0&0&0&-\frac{1}{2}&\frac{1}{2}&\frac{1}{2}&0\\\\[8.53581pt]
1&\frac{1}{4}&\frac{7}{4}&\frac{3}{4}&-\frac{7}{4}&0&0&\frac{3}{2}&-\frac{3}{2}&\frac{1}{2}&\frac{1}{2}\\\\[8.53581pt]
\hline\cr&\frac{1}{4}&\frac{7}{4}&\frac{3}{4}&-\frac{7}{4}&0&0&\frac{3}{2}&-\frac{3}{2}&\frac{1}{2}&\frac{1}{2}\end{array}.$

In space, a fifth order finite difference WENO scheme [29, 40], and a fourth
order compact central difference discretization, introduced in section 2.3,
are used to reconstruct the first and second order spatial derivatives,
respectively. Overall the resulting scheme is third order in time and fourth
order in space, and denoted as “T3S4”. For comparison, we also consider a
first order scheme, with a first order IMEX scheme (2.10) in time and a first
order numerical flux for (2.27), which are

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

while a second order central difference scheme [10] is used for second order
derivatives. The scheme is denoted as “T1S1”. Here the time step is taken to
be $\Delta t=\text{CFL}\,\Delta x/\Lambda$ with $\Lambda=\max\\{|{\bf
u}|+\min(1,1/\varepsilon)\sqrt{gh}\\}$ and $\text{CFL}=0.2$. Numerical results
will validate that our proposed scheme is high order in space and in time, AP
and AA in the diffusive limit, as well as well-balancing for capturing
propagation of small perturbations around a hydrostatic state.

######  Example 4.1.

(1D accuracy test with linear friction) In this example, we start to test the
order of convergence for the T3S4 scheme, with different values of
$\varepsilon$ in the one-dimensional case. For simplicity, we consider the
system (1.7) with a linear friction term where $\gamma=1$ is taken. We
consider a well-prepared initial condition, in order to show both AP and AA
properties, which is given as follows:

(4.2) $h(x,0)=\sin(\pi x)+2,\qquad q(x,0)=-2\pi\cos(\pi x)\left(\sin(\pi
x)+2\right),$

with a flat bottom $b(x)=0$. The computational domain is $[0,2]$ with a
periodic boundary condition, and $g=2$. We take three different values of
$\varepsilon$: $\varepsilon=1,10^{-2},10^{-6}$. Setting the final time
$T=0.02$, we compute the numerical errors and order of convergence, by
comparing to a reference solution with $N=2560$. The numerical errors and
orders are shown in Table 4.1. From the results, we can see that when
$\varepsilon=1$ and $\varepsilon=10^{-6}$, with mesh refinement, a 3rd order
accuracy has been obtained, namely the temporal accuracy dominates. For the
intermediate regime $\varepsilon=10^{-2}$, order reductions can be observed,
which is a typical case occurred in such high order IMEX methods [10, 11, 28].
From this table, we verify that the AP and AA properties are obtained.

Table 4.1.  Example 4.1. $L^{1}$ errors and orders for $h$ and $hu$ with different $\varepsilon$’s. $T=0.02$.  | N |  $1$ | $10^{-2}$ | $10^{-6}$
---|---|---|---|---
error | order | error | order | error | order
$h$ | 40 | 9.30E-04 | – | 4.88E-06 | – | 6.39E-06 | –
80 | 7.26E-05 | 3.68 | 2.73E-06 | 0.84 | 5.35E-07 | 3.58
160 | 2.18E-06 | 5.06 | 9.20E-07 | 1.57 | 5.89E-08 | 3.18
320 | 1.10E-07 | 4.31 | 2.39E-07 | 1.94 | 8.34E-09 | 2.82
640 | 7.02E-09 | 3.97 | 4.92E-08 | 2.28 | 1.11E-09 | 2.91
1280 | 5.62E-10 | 3.64 | 7.34E-09 | 2.75 | 1.26E-10 | 3.14
$hu$ | 40 | 6.35E-03 | – | 1.84E-04 | – | 1.02E-04 | –
80 | 4.64E-04 | 3.77 | 5.24E-05 | 1.81 | 8.49E-06 | 3.58
160 | 1.39E-05 | 5.06 | 1.48E-05 | 1.83 | 1.36E-06 | 2.65
320 | 6.78E-07 | 4.36 | 4.87E-06 | 1.60 | 2.02E-07 | 2.75
640 | 4.44E-08 | 3.93 | 1.20E-06 | 2.02 | 2.70E-08 | 2.90
1280 | 3.65E-09 | 3.61 | 2.01E-07 | 2.57 | 3.11E-09 | 3.11

######  Example 4.2.

(1D accuracy test with nonlinear friction) In this example, we will test the
order of convergence for our scheme with nonlinear friction. Following [52],
we construct an exact solution as follows:

(4.3) $h(x,t)=2+\varepsilon^{2}\sin(\pi(x-t)),\qquad
q(x,t)=2+\varepsilon^{2}\sin(\pi(x-t)).$

This solution satisfies the SWEs with extra source terms

(4.4) $\left\\{\begin{array}[]{ll}h_{t}+q_{x}=0,\\\\[8.53581pt] \displaystyle
q_{t}+\left(\frac{q^{2}}{h}\right)_{x}+\frac{1}{\varepsilon^{2}}\left(\frac{g}{2}h^{2}\right)_{x}=-\frac{1}{\varepsilon^{2}}gk^{2}\frac{|q|q}{h^{7/3}}+\frac{1}{\varepsilon^{2}}gk^{2}\left[2+\varepsilon^{2}\sin(\pi(x-t))\right]^{-1/3}\\\\[8.53581pt]
\hskip 125.19194pt+g\pi h\cos(\pi(x-t)),\end{array}\right.$

Here, a flat bottom $b(x)=0$ is considerred and the Manning coefficient is
$k=1$, and $g=1$. The computational domain is $\Omega=[0,2]$. For this
example, due to those extra source terms in (4.4), we are not able to take
$\varepsilon$ approaching $0$. Instead we set $\varepsilon=1$ and take a final
time $T=0.04$. We compute the numerical errors and orders by comparing to the
exact solutions, and show the results in Table 4.2. From the table, we observe
almost fifth-order accuracy. This might be due to that spatial errors of
convective terms are dominated in this hyperbolic regime.

Table 4.2.  Example 4.2. $L^{1}$ errors and orders for $h$ and $hu$ with $\varepsilon=1$. $T=0.04$. N |  $h$ | $hu$
---|---|---
error | order | error | order
10 | 2.37E-03 | – | 5.23E-03 | –
20 | 8.11E-05 | 4.87 | 1.51E-04 | 5.11
40 | 2.46E-06 | 5.04 | 3.04E-06 | 5.64
80 | 7.59E-08 | 5.02 | 7.89E-08 | 5.27
160 | 2.33E-09 | 5.02 | 2.05E-09 | 5.27
320 | 6.93E-11 | 5.07 | 4.13E-11 | 5.63
640 | 1.96E-12 | 5.14 | 4.72E-12 | 3.13
1280 | 1.71E-13 | 3.52 | 7.68E-13 | 2.62

######  Example 4.3.

In this example, we will verify the AP property of our scheme. We consider two
initial conditions: smooth and discontinuous, as in [13]. The spatial domain
for both cases is $[-5,5]$. The smooth initial condition is given by

(4.5) $h(x,0)=\left\\{\begin{aligned} &2,&x<-1;\\\
&\frac{1}{2}\left(3+\sin\left(\frac{3\pi x}{2}\right)\right),&-1\leq x<1;\\\
&1,&\text{otherwise};\end{aligned}\right.\qquad q(x,0)=0,$

while the discontinuous initial condition is:

(4.6) $h(x,0)=\left\\{\begin{aligned} &2,&x<0;\\\
&1,&\text{otherwise};\end{aligned}\right.\qquad q(x,0)=0.$

Since the stiffness of (1.7) is caused by the long time and strong friction,
to investigate this, we consider two tests: (1) we fix the final time $T=0.01$
and set the coefficient of friction $gk^{2}=1.0$, where $g=9.812$. We take
four different values of $\varepsilon$, $1,\,0.2,\,0.1,\,0.0005$, and compare
the numerical solutions with the solution obtained from the limiting equation
(2.7), which we denote as “lim”. The results obtained by T3S4 are shown in
Fig. 4.1, with $N=200$. From the figure, we observe that as $\varepsilon$ goes
to zero, the solutions approach to the solution for the limiting equation,
namely, our scheme preserve the AP property. (2) we fix $\varepsilon=1$
(namely without dimensionalization), but take a final time $T=0.01\theta$
where $\theta=\sqrt{gk^{2}}$, to test the long-time behavior of the solutions.
We take four different values of $\theta=1,\,5,\,10,\,50$. The results of
water depth $h$ are shown in Fig. 4.2 with the same grid $N=200$. It can be
seen that the numerical results are almost the same as the first case, which
confirm the scaling in (1.2).

Finally, we compare the results obtained by the high order T3S4 scheme with
those obtained by the first order T1S1 scheme. We take two cases. In the first
case, the coefficient of the friction is $gk^{2}=1$ and the final time is
$T=0.01$, where the parameter $\varepsilon=0.0005$. In the second case, we
choose $gk^{2}=2500$, $T=0.5$, and $\varepsilon=1$. The results are shown in
Fig. 4.3 and Fig. 4.4, respectively. From the results, we observe that as
$\varepsilon$ goes to zero, or the simulation time increases, the first order
solution only match the high order solution under a very refined mesh size,
which show the necessity of using high order schemes.

(a) water depth $h$

(b) momentum $q$

(c) water depth $h$

(d) momentum $q$

Figure 4.1.  Example 4.3. Numerical solutions of the water depth $h$ (left)
and the momentum $q$ (right) for the smooth initial condition (4.5) (top) and
the discontinuous initial condition (4.6) (bottom) with different
$\varepsilon$’s .

(a) water depth $h$

(b) water depth $h$

Figure 4.2.  Example 4.3. Numerical solutions of the water depth $h$ for the
smooth initial condition (4.5) (left) and the discontinuous initial condition
(4.6) (right) at different time.

(a) $\varepsilon=0.0005$, $gk^{2}=1$, $T=0.01$

(b) $\varepsilon=0.0005$, $gk^{2}=1$, $T=0.01$

Figure 4.3.  Example 4.3. Numerical solutions of the water depth $h$ for the
smooth initial condition (4.5) (left) and the discontinuous initial condition
(4.6) (right), with the T1S1 and T3S4 schemes.

(a) $\varepsilon=1$, $gk^{2}=2500$, $T=0.5$

(b) $\varepsilon=1$, $gk^{2}=2500$, $T=0.5$

Figure 4.4.  Example 4.3. Numerical solutions of the water depth $h$ for the
smooth initial conditions (4.5) (left) and the discontinuous initial condition
(4.6) (right), with the T1S1 and T3S4 schemes.

######  Example 4.4.

In this example, we now test the ability of our proposed scheme to capture the
propagation of small perturbations over strong friction. The initial condition
is taken to be:

(4.7) $h(x,0)=\left\\{\begin{aligned} &0.5+0.1\sin(10\pi(x-1.1)),&1.1\leq
x\leq 1.2;\\\ &0.5,&\text{otherwise};\end{aligned}\right.\quad q(x,0)=0.$

The computational domain is $[0,2]$, with a gravitational constant $g=9.812$
and a flat bottom $b(x)=0$. We set the friction coefficient be $k=1$,
corresponding to strong friction. Inflow and outflow boundary conditions are
considerred. We compute the solution up to a final time $T=0.02$. Five
different values of $\varepsilon$ are taken, i.e
$1,\,0.5,\,0.1,\,10^{-2},\,10^{-6}$. As $\varepsilon$ becoming smaller, the
friction becomes stronger. We compare the results with the one obtained by the
limiting equation (2.7), which is denoted as “lim”, and present the numerical
results in Fig. 4.5. From the results, we observe that as $\varepsilon$
decreases, the diffusion of the system becomes to dominate. Our scheme can be
uniformly stable and capture the small perturbations well. Besides, when
$\varepsilon=10^{-6}$, the solution matches very well with the one obtained by
the limiting equation (2.7), namely our scheme is AP in the diffusive regime.

(a) the surface level $H$

(b) momentum $q$

Figure 4.5.  Example 4.4. Numerical solutions of the surface level $H$ (left)
and the momentum $q$ (right), with strong friction.

######  Example 4.5.

In this example, we will check the well-balanced property of our scheme. An
equilibrium initial condition is taken, with a small perturbation adding to a
surface level, following the study in [52]. This setup is an extension of the
work presented in [46, 49, 28], where the system is considered without
friction.

Here we consider a non-flat trigonometric bottom topology, which is given by

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

The initial condition with a small disturbance is taken:

(4.9) $h(x,0)+b(x)=\left\\{\begin{aligned} &1.001,&1.1\leq x\leq 1.2;\\\
&1,&\text{otherwise};\end{aligned}\right.\qquad q(x,0)=0.$

Here, the computational domain is $[0,2]$, $g=9.812$ and $k=1$. The final time
$T=0.2$.

To demonstrate that our scheme is not affected by the dimensionless parameter
$\varepsilon$, we scale the initial condition into a dimensionless form, using
$\hat{t}=\varepsilon t$, $\hat{k}=\varepsilon k$, and
$\hat{q}=\frac{q}{\varepsilon}$, as shown in (1.2) and (1.5). We take $N=200$,
and consider an inflow and outflow boundary condition. Two different values of
$\varepsilon$ are used, i.e., $1$ and $0.6$. After we obtain the numerical
results for different $\varepsilon$’s, we transform them back to the dimension
case. The results are shown in Fig. 4.6. We can see the results match each
other very well, namely, they are not affected by the scaling. Besides, these
results are also the same as in [52], which can capture the small
perturbations well. This demonstrates that our scheme is well-balanced
independent of the parameter $\varepsilon$. Finally, we compare the results
obtained by the high order T3S4 scheme and the first order T1S1 scheme, which
are shown in Fig. 4.7. Clearly the high order results are much better than the
first order results, and they match well when the mesh is fine enough for the
first order scheme.

Figure 4.6.  Example 4.5. Numerical solutions of the surface level $H$ (left)
and the momentum $q$ (right) for the T3S4 scheme with different
$\varepsilon$’s.

Figure 4.7.  Example 4.5. Numerical solutions of the surface level $H$ (left)
and the momentum $q$ (right), with the T1S1 and T3S4 schemes. $\varepsilon=1$.

######  Example 4.6.

In this example, we consider a shock wave with a non-flat bottom [52]. The
Manning coefficient is $k=0.09$, $g=9.812$. The non-flat bottom has two humps
which is defined as follows:

(4.10) $b(x)=\left\\{\begin{aligned} &\max(0,0.25-20(x-0.25)^{2}),&0\leq x\leq
0.5;\\\ &\max(0,0.25-20(x-0.75)^{2}),&0.5<x\leq 1;\end{aligned}\right.$

The initial condition is given by

(4.11) $h(x,0)+b(x)=\left\\{\begin{aligned} &1,&0\leq x\leq 0.5;\\\
&0.5,&0.5<x\leq 1;\end{aligned}\right.\quad q(x,0)=0.$

We take a computational domain $[0,1]$ with $N=200$ and set $\varepsilon=1$.
An inflow and outflow boundary condition is adopted, with a final time
$T=0.08$. The numerical results are shown in Fig. 4.8. From the results, we
observe that two waves are generated from the initial discontinuity and
propagate to left and right, respectively. We also compare the results
obtained by the high order T3S4 scheme, with those obtained by the first order
T1S1 scheme. Similarly the high order scheme has better resolution than the
first order scheme. This example demonstrates the good performance of our high
order scheme for shock capturing.

(a) water surface $h+b$

(b) momentum $q$

Figure 4.8.  Example 4.6. Numerical solutions of the surface level $H$ (left)
and the momentum $q$ (right). $\varepsilon=1$. $T=0.08$.

######  Example 4.7.

(2D accuracy test with linear friction) In this example, we test the orders of
accuracy for our scheme in the 2D case. Similarly we consider the system (1.7)
with a linear friction of $\gamma=1$. The initial conditions are given by:

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

The bottom is flat with $b(x,y)=0$. The computational domain is
$\Omega=[0,2]^{2}$, and we set $g=2$. Periodic boundary conditions are used.
The computational domain is divided uniformly into $N_{\ell}^{2}$ meshes,
where $N_{\ell}=N_{0}\cdot 2^{k}$ with $N_{0}=8$ and $\ell=0,1,2,3,4,5$. Since
the exact solution is not available, we take the numerical solution obtained
by $N^{2}=512^{2}$ as the reference solution, and compute the numerical errors
by comparing the numerical solution and the reference solution. Here we
consider three different values of $\varepsilon$: $1,10^{-2},10^{-6}$.

The numerical errors and orders are shown in Table 4.3 at a final time
$T=0.01$. From the results, we observe that around fourth order accuracy for
$\varepsilon=1$, and third order for $\varepsilon=10^{-6}$ are obtained, which
is similar to the one-dimensional case. This demonstrates that our scheme
performs well in both hyperbolic and diffusive regimes. An order reduction
also appears for the intermediate regime of $\varepsilon=10^{-2}$, similar to
the accuracy tests in [10, 11, 28].

Table 4.3.  Example 4.7. $L^{1}$ errors and orders for $h$, $hu$, and $hv$ with $\varepsilon=1,10^{-2},10^{-6}$. | N |  $1$ | $10^{-2}$ | $10^{-6}$
---|---|---|---|---
error | order | error | order | error | order
$h$ | 16 | 5.94E-03 | – | 2.57E-04 | – | 2.97E-04 | –
32 | 5.62E-04 | 3.40 | 1.25E-05 | 4.36 | 1.78E-05 | 4.06
64 | 2.37E-05 | 4.57 | 4.50E-06 | 1.48 | 1.37E-06 | 3.71
128 | 1.31E-06 | 4.17 | 1.40E-06 | 1.69 | 1.32E-07 | 3.37
256 | 8.14E-08 | 4.01 | 2.78E-07 | 2.33 | 1.39E-08 | 3.25
$hu$ | 16 | 3.82E-02 | – | 5.16E-03 | – | 5.52E-03 | –
32 | 3.72E-03 | 3.36 | 3.14E-04 | 4.04 | 2.94E-04 | 4.23
64 | 1.67E-04 | 4.48 | 9.44E-05 | 1.73 | 1.80E-05 | 4.03
128 | 9.43E-06 | 4.15 | 3.00E-05 | 1.65 | 2.47E-06 | 2.86
256 | 5.52E-07 | 4.10 | 7.29E-06 | 2.04 | 3.31E-07 | 2.90
$hv$ | 16 | 3.82E-02 | – | 5.16E-03 | – | 5.52E-03 | –
32 | 3.72E-03 | 3.36 | 3.14E-04 | 4.04 | 2.94E-04 | 4.23
64 | 1.67E-04 | 4.48 | 9.44E-05 | 1.73 | 1.80E-05 | 4.03
128 | 9.43E-06 | 4.15 | 3.00E-05 | 1.65 | 2.47E-06 | 2.86
256 | 5.52E-07 | 4.10 | 7.29E-06 | 2.04 | 3.31E-07 | 2.90

######  Example 4.8.

In this example, we now test the accuracy of the scheme for nonlinear friction
in the two-dimensional case. Similar to the one-dimensional problem, the
“manufactured” exact solutions are taken to be:

$\displaystyle h(x,y,t)=2+\varepsilon^{2}\sin(\pi(x+y-2t)),$ $\displaystyle
q_{1}(x,y,t)=2+\varepsilon^{2}\sin(\pi(x+y-2t)),$ $\displaystyle
q_{2}(x,y,t)=2+\varepsilon^{2}\sin(\pi(x+y-2t)),$

which satisfy the SWEs with extra source terms

(4.14)
$\left\\{\begin{array}[]{ll}\partial_{t}h+\partial_{x}q_{1}+\partial_{y}q_{2}=0,\\\\[8.53581pt]
\partial_{t}q_{1}+\partial_{x}\left(\frac{q_{1}^{2}}{h}\right)+\partial_{y}\left(\frac{q_{1}q_{2}}{h}\right)+\frac{1}{\varepsilon^{2}}\partial_{x}\left(\frac{g}{2}h^{2}\right)=-\frac{1}{\varepsilon^{2}}\frac{gk^{2}\sqrt{q_{1}^{2}+q_{2}^{2}}}{h^{7/3}}q_{1}+\frac{gk^{2}}{\varepsilon^{2}}\left[2+\right.\\\\[8.53581pt]
\hskip 91.04872pt\left.+\varepsilon^{2}\sin(\pi(x+y-2t))\right]^{-1/3}+g\pi
h\cos(\pi(x+y-2t)),\\\\[8.53581pt]
\partial_{t}q_{2}+\partial_{x}\left(\frac{q_{1}q_{2}}{h}\right)+\partial_{y}\left(\frac{q_{2}^{2}}{h}\right)+\frac{1}{\varepsilon^{2}}\partial_{y}\left(\frac{g}{2}h^{2}\right)=-\frac{1}{\varepsilon^{2}}\frac{gk^{2}\sqrt{q_{1}^{2}+q_{2}^{2}}}{h^{7/3}}q_{2}+\frac{gk^{2}}{\varepsilon^{2}}\left[2+\right.\\\\[8.53581pt]
\hskip 91.04872pt\left.+\varepsilon^{2}\sin(\pi(x+y-2t))\right]^{-1/3}+g\pi
h\cos(\pi(x+y-2t)).\end{array}\right.$

$g=1$ and $k=1$. The computational domain is $[0,2]^{2}$, and the final time
is $T=0.04$. We only consider $\varepsilon=1$ as the one-dimensional case. The
errors and orders are shown in Table 4.4. From the table, we observe a desired
fifth order accuracy.

Table 4.4.  Example 4.8. The $L^{1}$ errors and orders for $h$, $hu$, and $hv$ with $\varepsilon=1$. N |  $h$ | $hu$ | $hv$
---|---|---|---
error | order | error | order | error | order
8 | 3.20E-02 | – | 6.46E-02 | – | 6.46E-02 | –
16 | 1.74E-03 | 4.20 | 3.62E-03 | 4.16 | 3.62E-03 | 4.16
32 | 6.36E-05 | 4.78 | 1.03E-04 | 5.13 | 1.03E-04 | 5.13
64 | 2.00E-06 | 4.99 | 3.12E-06 | 5.05 | 3.12E-06 | 5.05
128 | 6.23E-08 | 5.01 | 9.62E-08 | 5.02 | 9.62E-08 | 5.02
256 | 1.93E-09 | 5.01 | 2.93E-09 | 5.04 | 2.93E-09 | 5.04

######  Example 4.9.

In the final example, we check the well-balanced property for our scheme in
the two-dimensional case. We consider a bottom topography being an isolated
elliptical shaped hump

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

and the initial conditions are

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

This example has been studied in [52]. The computational domain is
$[0,2]\times[0,1]$, with an inflow and outflow boundary condition in the $x$
direction and a periodic boundary condition in the $y$ direction. Similar
experiments were conducted in prior works such as [46, 49, 28], without
friction. In this case, a small friction has been added with Manning
coefficient $k=0.09$. The numerical solutions of the surface level $H=h+b$ are
shown in Fig. 4.9 with two different meshes, i.e., $200\times 100$ and
$400\times 200$. The initial perturbation is separated into two waves
propagating to 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 our proposed T3S4 scheme. The numerical
results are similar to those in [52].

(a) $t=0.12$

(b) $t=0.12$

(c) $t=0.24$

(d) $t=0.24$

(e) $t=0.36$

(f) $t=0.36$

(g) $t=0.48$

(h) $t=0.48$

(i) $t=0.6$

(j) $t=0.6$

Figure 4.9.  Example 4.9. Numerical solutions of the surface level $H$. From
top to bottom: at $t=0.12$ from 0.9998 to 1.0060; at $t=0.24$ from 0.9967 to
1.0130; at $t=0.36$ from 0.9901 to 1.0097; $t=0.48$ from 0.9906 to 1.0043;
$t=0.6$ from 0.9955 to 1.0045, 30 contour lines are used. Mesh: $200\times
100$ (left); $400\times 200$ (right).

##  5\. Conclusion

In this paper, we have proposed a high order AP finite difference WENO scheme
for the SWEs with Manning friction and a bottom topology. The scheme is proved
to be AP, AA and well-balanced, which are also confirmed by one- and two-
dimensional numerical tests. In this work, the diffusion limit are added as a
penalization, and an additive IMEX RK method has been utilized. In the future,
we may further consider a semi-implicit approach without adding the diffusion
terms.

## Acknowledgements

Guanlan Huang and Tao Xiong are partially supported by National Key R&D
Program of China No. 2022YFA1004500, NSFC grant No. 11971025 and 92270112, NSF
of Fujian Province No. 2023J02003, and the Strategic Priority Research Program
of Chinese Academy of Sciences Grant No. XDA25010401.

Sebastiano Boscarino is partially supported by the Italian Ministry of
Instruction, University and Research (MIUR), with funds coming from PRIN
Project 2022 (2022KA3JBA, entitled “Advanced numerical methods for time
dependent parametric partial differential equations and applications”) and
from Italian Ministerial grant PRIN 2022 PNRR “FIN4GEO: Forward and Inverse
Numerical Modeling of hydrothermal systems in volcanic regions with
application to geothermal energy exploitation.”, No. P2022BNB97. Furthermore,
Sebastiano Boscarino is supported for this work by the Spoke 1 ”FutureHPC &
BigData” of the Italian Research Center on High-Performance Computing, Big
Data and Quantum Computing (ICSC) funded by MUR ”Missione 4 Componente 2
Investimento 1.4: Potenziamento strutture di ricerca e creazione di campioni
nazionali di R&S (M4C2-19 )”. S. Boscarino is a member of the INdAM Research
group GNCS.

## References

  * [1] F. Andreu, J. Mazon, S. de León, and J. Toledo.  Existence and uniqueness for a degenerate parabolic equation with ${L}^{1}$-data.  Transactions of the American Mathematical Society, 351(1):285–306, 1999.
  * [2] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri.  Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations.  Applied Numerical Mathematics, 25(2):151–167, 1997.
  * [3] J. W. Barrett and W. Liu.  Finite element approximation of the parabolic p-Laplacian.  SIAM Journal on Numerical Analysis, 31(2):413–428, 1994.
  * [4] J. W. Barrett and W. B. Liu.  Finite element approximation of the p-Laplacian.  Mathematics of Computation, 61(204):523–537, 1993.
  * [5] A. Bermudez and M. E. Vazquez.  Upwind methods for hyperbolic conservation laws with source terms.  Computers & Fluids, 23(8):1049–1071, 1994.
  * [6] C. Berthon, P. LeFloch, and R. Turpault.  Late-time/stiff-relaxation asymptotic-preserving approximations of hyperbolic equations.  Mathematics of Computation, 82(282):831–860, 2013.
  * [7] S. Boscarino, F. Filbet, and G. Russo.  High order semi-implicit schemes for time dependent partial differential equations.  Journal of Scientific Computing, 68(3):975–1001, 2016.
  * [8] S. Boscarino, P. G. LeFloch, and G. Russo.  High-order asymptotic-preserving methods for fully nonlinear relaxation problems.  SIAM Journal on Scientific Computing, 36(2):A377–A395, 2014.
  * [9] 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.
  * [10] 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.
  * [11] 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.  SIAM Journal on Scientific Computing, 4:B368–B394, 2022.
  * [12] S. Boscarino and G. Russo.  Flux-explicit IMEX Runge-Kutta schemes for hyperbolic to parabolic relaxation problems.  SIAM Journal on Numerical Analysis, 51(1):163–190, 2013.
  * [13] S. Bulteau, M. Badsi, C. Berthon, and M. Bessemoulin-Chatard.  A fully well-balanced and asymptotic preserving scheme for the shallow-water equations with a generalized Manning friction source term.  Calcolo, 58(4):41, 2021.
  * [14] J. C. Butcher.  Numerical Methods for Ordinary Differential Equations.  Third Edition, John Wiley & Sons Ltd, 2016.
  * [15] 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.
  * [16] G.-Q. Chen, C. D. Levermore, and T.-P. Liu.  Hyperbolic conservation laws with stiff relaxation terms and entropy.  Communications on Pure and Applied Mathematics, 47(6):787–830, 1994\.
  * [17] A. Chertock, S. Cui, A. Kurganov, and T. Wu.  Well-balanced positivity preserving central-upwind scheme for the shallow water system with friction terms.  International Journal for Numerical Methods in Fluids, 78(6):355–383, 2015.
  * [18] F. Crispo, P. Maremonti, and M. Růžička.  Global $L^{r}$-estimates and regularizing effect for solutions to the $p(t,x)$-Laplacian systems.  Advances in Differential Equations, 24(7/8):407 – 434, 2019.
  * [19] S. S. de León and J. Toledo.  Regularity for entropy solutions of parabolic p-Laplacian type equations.  Publicacions Matematiques, pages 665–683, 1999.
  * [20] J. I. Díaz and M. Herrero.  Estimates on the support of the solutions of some nonlinear elliptic and parabolic problems.  Proceedings of the Royal Society of Edinburgh Section A: Mathematics, 89(3-4):249–258, 1981.
  * [21] J. Eom and R. Sato.  Large time behavior of ode type solutions to parabolic p-laplacian type equations.  Communications on Pure & Applied Analysis, 19(9), 2020.
  * [22] A. Ern, S. Piperno, and K. Djadel.  A well-balanced Runge–Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying.  International Journal for Numerical Methods in Fluids, 58(1):1–25, 2008.
  * [23] J. García Azorero and I. Peral Alonso.  Existence and nonuniqueness for the p-Laplacian.  Communications in Partial Differential Equations, 12(12):126–202, 1987.
  * [24] I. Gómez-Bueno, M. J. C. Díaz, C. Parés, and G. Russo.  Collocation methods for high-order well-balanced methods for systems of balance laws.  Mathematics, 9(15):1799, 2021.
  * [25] E. Hairer and G. Wanner.  Solving ordinary differential equations II: stiff and differential algebraic problems, volume 2.  Springer Verlag, 1993.
  * [26] M. A. Herrero and J. L. Vázquez.  Asymptotic behaviour of the solutions of a strongly nonlinear parabolic problem.  In Annales de la Faculté des sciences de Toulouse: Mathématiques, volume 3, pages 113–127, 1981.
  * [27] J. Hu, S. Jin, and Q. Li.  Asymptotic-preserving schemes for multiscale hyperbolic and kinetic equations.  Handbook of Numerical Analysis, 18:103–129, 2017.
  * [28] G. Huang, Y. Xing, and T. Xiong.  High order well-balanced asymptotic preserving finite difference WENO schemes for the shallow water equations in all Froude numbers.  Journal of Computational Physics, 463:111255, 2022.
  * [29] G.-S. Jiang and C.-W. Shu.  Efficient implementation of weighted ENO schemes.  Journal of Computational Physics, 126(1):202–228, 1996.
  * [30] S. Jin.  Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations.  SIAM Journal on Scientific Computing, 21(2):441–454, 1999.
  * [31] S. Jin.  Asymptotic-preserving schemes for multiscale physical problems.  Acta Numerica, 31:415–489, 2022.
  * [32] S. Kamin and J. L. Vázquez.  Fundamental solutions and asymptotic behaviour for the $p$-Laplacian equation.  Revista Matemática Iberoamericana, 4(2):339–354, 1988.
  * [33] A. Kurganov.  Finite-volume schemes for shallow-water equations.  Acta Numerica, 27:289–351, 2018.
  * [34] A. Kurganov and D. Levy.  Central-upwind schemes for the Saint-Venant system.  ESAIM: Mathematical Modelling and Numerical Analysis, 36(3):397–425, 2002.
  * [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] V. Michel-Dansac, C. Berthon, S. Clain, and F. Foucher.  A well-balanced scheme for the shallow-water equations with topography or Manning friction.  Journal of Computational Physics, 335:115–154, 2017.
  * [37] 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.
  * [38] Z. Peng, Y. Cheng, J.-M. Qiu, and F. Li.  Stability-enhanced AP IMEX-LDG schemes for linear kinetic transport equations under a diffusive scaling.  Journal of Computational Physics, 415:109485, 2020.
  * [39] S. Rhebergen, O. Bokhove, and J. J. van der Vegt.  Discontinuous galerkin finite element methods for hyperbolic nonconservative partial differential equations.  Journal of Computational Physics, 227(3):1887–1922, 2008.
  * [40] C.-W. Shu.  Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws.  Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, pages 325–432, 1998.
  * [41] L. Song, J. Zhou, Q. Li, X. Yang, and Y. Zhang.  An unstructured finite volume model for dam-break floods with wet/dry fronts over complex topography.  International Journal for Numerical Methods in Fluids, 67(8):960–980, 2011.
  * [42] K. Wu and Y. Xing.  Uniformly high-order structure-preserving discontinuous Galerkin methods for Euler equations with gravitation: Positivity and well-balancedness.  SIAM Journal on Scientific Computing, 43(1):A472–A510, 2021.
  * [43] X. Xia and Q. Liang.  A new efficient implicit scheme for discretising the stiff friction terms in the shallow water equations.  Advances in Water Resources, 117:87–97, 2018.
  * [44] X. Xia, Q. Liang, X. Ming, and J. Hou.  An efficient and stable hydrodynamic model with novel source term discretization schemes for overland flow and flood simulations.  Water Resources Research, 53(5):3730–3759, 2017.
  * [45] Y. Xing.  Numerical methods for the nonlinear shallow water equations.  In R. Abgrall and C.-W. Shu, editors, Handbook of Numerical Analysis, Volume 18, Handbook of Numerical Methods for Hyperbolic Problems: Applied and Modern Issues, volume 18, pages 361–384. North-Holland, Elsevier, Amsterdam, 2017.
  * [46] 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.
  * [47] Y. Xing and C.-W. Shu.  A new approach of high order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms.  Communications in Computational Physics, 1(1):100–134, 2006.
  * [48] Y. Xing and C.-W. Shu.  High-order well-balanced finite difference weno schemes for a class of hyperbolic systems with source terms.  Journal of Scientific Computing, 27:477–494, 2006.
  * [49] Y. Xing and C.-W. Shu.  High order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms.  Journal of Computational Physics, 214(2):567–598, 2006.
  * [50] 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.
  * [51] T. Xiong, W. Sun, Y. Shi, and P. Song.  High order asymptotic preserving discontinuous Galerkin methods for gray radiative transfer equations.  Journal of Computational Physics, 463:111308, 2022.
  * [52] R. Yang, Y. Yang, and Y. Xing.  High order sign-preserving and well-balanced exponential Runge-Kutta discontinuous Galerkin methods for the shallow water equations with friction.  Journal of Computational Physics, 444:110543, 2021.
  * [53] J. G. Zhou, D. M. Causon, C. G. Mingham, and D. M. Ingram.  The surface gradient method for the treatment of source terms in the shallow-water equations.  Journal of Computational Physics, 168(1):1–25, 2001.