# Adaptive Central-Upwind Scheme on Triangular Grids for the Shallow Water
Model with variable density

Thuong Nguyen111Department of Mathematics, The University of Utah, Salt Lake
City, UT 84112, USA<EMAIL_ADDRESS>
###### Abstract

In this paper, we construct a robust adaptive central-upwind scheme on
unstructured triangular grids for two-dimensional shallow water equations with
variable density. The method is well-balanced, positivity-preserving, and
oscillation free at the curve where two types of fluid merge. The proposed
approach is an extension of the adaptive well-balanced, positivity-preserving
scheme developed in Epshteyn and Nguyen (arXiv preprint arXiv:2011.06143,
2020). In particular, to preserve “lake-at-rest” steady states, we utilize the
Riemann Solver with appropriately rotated coordinates to obtain the point
values in neighborhood of the fluid interface. In addition, to improve the
efficiency of an adaptive method in the multifluid flow, the curve of density
discontinuity is reconstructed by using the level set method and volume
fraction method. To demonstrate the accuracy, high-resolution, and efficiency
of the new adaptive central-upwind scheme, several challenging tests for
Shallow water models with variable density are performed.

###### keywords:

Shallow water equations with variable density, central-upwind scheme, well-
balanced and positivity-preserving scheme, adaptive algorithm, interface
tracking, Riemann solver, weak local residual error estimator, unstructured
triangular grid

AMS subject classification: 76M12, 65M08, 35L65, 86-08, 86A05

##  1 Introduction

The main goal of this paper is to develop an adaptive well-balanced
positivity-preserving central-upwind scheme on triangular grids for shallow
water equations with variable density (SWEDs). The two-dimensional (2-D)
system of SWEDs can be written as,

$\displaystyle w_{t}+(hu)_{x}+(hv)_{y}=0,$ (1.1a)
$\displaystyle(hu)_{t}+\Big{(}hu^{2}+\frac{g}{2\rho_{0}}h^{2}\rho\Big{)}_{x}+(huv)_{y}=-\frac{g}{\rho_{0}}h\rho
B_{x},$ (1.1b)
$\displaystyle(hv)_{t}+(huv)_{x}+\Big{(}hv^{2}+\frac{g}{2\rho_{0}}h^{2}\rho\Big{)}_{y}=-\frac{g}{\rho_{0}}h\rho
B_{y},$ (1.1c) $\displaystyle(h\rho)_{t}+(hu\rho)_{x}+(hv\rho)_{y}=0,$ (1.1d)

where $t$ is the time, $x$ and $y$ are spatial coordinates ($(x,y)\in\Omega$),
$h(x,y,t)$ is the water height, $\rho(x,y,t)$ is the density, $u(x,y,t)$ and
$v(x,y,t)$ are the $x$\- and $y$-components of the flow velocity, $B(x,y)$ is
the bottom topography, $g$ is the constant gravitational acceleration, and
$\rho_{0}$ is the reference density. The system (1.1a)–(1.1d) was proposed in
[10, 35, 36, 15] as a variation of the Saint-Venant equations to model multi-
phase flows in estuaries or deep ocean currents. The derivation of the system
is based on hydrostatic approximation which eliminates the variability in the
z-direction. The design of robust and accurate numerical algorithms for
computing the solutions of SWEDs system is an important and challenging
problem that has been extensively studied in the recent years.

A number of numerical schemes for balance laws have been introduced in recent
years, [28, 24, 20, 5, 27, 7, 6, 8, 23, 22, 2, 26, 3, 4, 39, 29]. Most of them
utilize a Riemann problem solver for the upwind evolution of the calculated
solution. However, as discussed in [9], the eigensystem of the system
(1.1a)–(1.1d) may be incomplete due to the resonance phenomenon. Hence, it may
be very difficult to design a reliable upwind scheme for the SWEDs. In our
paper, we therefore use central-upwind schemes which are Riemann-problem-
solver free methods, [21, 25, 29]. Central-upwind schemes have been referred
to “black-box” solvers for general multidimensional systems of hyperbolic
systems of conservation laws. In our prior work [13], we have derived a
successful adaptive central-upwind method for Saint-Venant system on
triangular grids.

Similar to the Saint-Venant system, a good method for SWEDs system should
preserve the non-negativity of $h$ and $\rho$, which is called the positivity-
preserving property. In addition, the scheme must ensure a well-balanced
property obtained when the numerical method preserve “lake-at-rest” steady-
state solutions. Otherwise, the numerical method may lead to significant
oscillations. Note that, the system (1.1a)–(1.1d) admits the following two
“lake-at-rest” steady-state solutions, [9]:

$w=\max\big{\\{}C,B(x,y)\big{\\}},\quad C=\rm{Const},\quad\rho=P\equiv
Const,\quad u\equiv v\equiv 0,$ (1.2)

and

$B\equiv\rm{Const},\quad h^{2}\rho\equiv Const,\quad u\equiv v\equiv 0.$ (1.3)

Preserving the solution (1.3) is a big challenge for numerically solving the
system (1.1a)–(1.1d) since using the conventional central-upwind methods may
not ensure the variable $h^{2}\rho$, so-called variable pressure, to be
constant at the contact waves. Therefore, for our adaptive scheme of SWEDs, we
will utilize the central-upwind scheme which is derived from [9] for shallow
water model with horizontal temperature variable (the system in [9] has
similar properties with the SWEDs). In [9], the proposed second-order semi-
discrete central-upwind scheme is capable of preserving the “lake at rest”
steady state (1.2) and (1.3) as well as the positivity of the water depth and
the temperature (the temperature variable is equivalent to the variable
density in our work). In particular, to preserve the second type of “lake at
rest” steady state (1.3) and suppress the pressure oscillations across the
interface, an efficient interface tracking method is performed. The main idea
of the interface approach in [9] is to completely avoid to use the information
from the cells where two types of fluids are numerically mixed, so-called
“mixed” cells when evolving the solution in the neighborhood of the interface.
The data in the “mixed” cells is replaced by the interpolated values that are
calculated using the reliable information from the nearby “single fluid”
cells. Namely, the point values in “mixed” cells are obtained by using the
approximated solution of the 1-D Riemann problems between the reliable “single
fluid” cell averages. However, the central-upwind method and the interface
tracking in [9] are designed on structured rectangular grids. In practice, one
needs to deal with complicated geometries, where the use of triangular grids
could be advantageous or even unavoidable. Hence, in this study, we extend the
interface tracking method from [9] to unstructured triangular grids by
utilizing the idea of the rotated coordinates proposed in [1]. The 1-D Riemann
solver [9] is then performed in the direction of normal vectors on each edge
of “mixed” cells. This extended central-upwind scheme maintains the well-
balanced and positivity-preserving properties.

In addition to achieving the accuracy of the solution, minimizing
computational cost is also one of the major challenge in the modeling and
numerical analysis of hydrodynamics. The traditional numerical schemes for the
system (1.1a)–(1.1d) are based on the use of very fine fixed meshes to
reconstruct delicate features of the solution. However, this can lead to high
computational cost, as well as poor resolution of all small scale features of
the problem. In many engineering and scientific applications, it is beneficial
to use adaptive meshes for improving the accuracy of the approximation at a
much lower cost. Therefore, another goal of this work is to design an adaptive
numerical algorithm for shallow water equations with variable density. There
is some very recent effort on the design of adaptive well-balanced and
positivity-preserving central-upwind schemes on quad-tree grids for shallow
water models [37, 15, 14], but no research has been done for the development
of such adaptive schemes for the shallow water with variable density on
unstructured triangular grids.

As a part of the mesh reconstruction, we will need to project the data from
the old mesh onto a new adaptive mesh. Since we avoid using the cell averages
of mixed cells, the data of a new cell which is reconstructed from an old
mixed cell is calculated based on the location of the density jumps and the
data from the nearby “reliable” cells. Therefore, the interface separating
different fluid phases must be accurately tracked or captured. In each “mixed”
cell, we will reconstruct an approximate interface, which is called “interface
reconstruction”. Many methods have been developed for this purpose in the last
two decades, [31, 38, 17, 34, 42, 41]. In our work, we consider the interface
reconstruction method derived in [45] due to the great advantages of this
approach such as high-resolution, efficiency and simplicity.

This paper is organized as follows. In Section 2, we present a well-balanced
positivity-preserving central-upwind scheme on unstructured triangular grids
for SWEDs which serves as the underlying discretization for the developed
adaptive algorithm. Next, in Section 3.1, we present the procedure to detect
mixed cells. The interface approximation will be briefly reviewed in Section
3.2. In Section 3.3, we discuss the cell averages correction, which is used to
prevent the density diffusion around the density jumps. We summarize the
adaptive central-upwind method in Section 4.1. We discuss a strategy of
adaptive mesh refinement in Section 4.2. In Section 4.3, we present an
adaptive second-order strong stability preserving Runge-Kutta method, employed
as a part of time evolution for the adaptive central-upwind scheme. We develop
a local posteriori error estimator in Section 4.4 which is used as a robust
indicator for the adaptive mesh refinement in our work. Finally, in Section 5,
we illustrate the high accuracy and efficiency of the developed adaptive
central-upwind scheme on a number of challenging tests for multi-fluid shallow
water models.

##  2 The Central Upwind Scheme for SWEDs in Triangular Mesh

In this section, we focus on developing the central-upwind method for the
SWEDs system (1.1a)-(1.1d), [15]. In particular, we will extend the adaptive
scheme from the Saint-Venant system in [13] to the SWEDs system. The developed
scheme will eliminate the oscillation appearing at the interface and ensure
the well-balanced and positivity-preserving properties.

We first rewrite the SWEDs system (1.1a)-(1.1d) into the vector form as,

$\mathbf{U}_{t}+\mathbf{F}(\mathbf{U},B)_{x}+\mathbf{G}(\mathbf{U},B)_{y}=\mathbf{S}(\mathbf{q},B),$
(2.1)

where

$U=(w,hu,hv,h\rho),$

and the fluxes and source term are:

$\displaystyle\mathbf{F}(\mathbf{U},B)$
$\displaystyle=(hu,\dfrac{(hu)^{2}}{w-B}+\dfrac{g}{2\rho_{0}}\rho(w-B)^{2},\dfrac{(hu)(hv)}{w-B},hu\rho)^{T},$
(2.2) $\displaystyle\mathbf{G}(\mathbf{q},B)$
$\displaystyle=(hv,\dfrac{(hu)(hv)}{w-B},\dfrac{(hv)^{2}}{w-B}+\dfrac{g}{2\rho_{0}}\rho(w-B)^{2},hv\rho)^{T},$
$\displaystyle\mathbf{S}(\mathbf{U},B)$
$\displaystyle=(0,-g\frac{\rho}{\rho_{0}}(w-B)B_{x},-g\frac{\rho}{\rho_{0}}(w-B)B_{y},0)^{T}.$

In our research, we consider the triangular mesh illustrated in Fig. 2.1 with
the following notations.

Figure 2.1:  A typical triangular cell with three neighbors.

${\mathcal{T}}:=\bigcup_{j}T_{j}$ is an unstructured triangulation of the
computational domain $\Omega$;

$T_{j}\in\mathcal{T}$ is a triangular cell of size $|T_{j}|$ with the
barycenter $(x_{j},y_{j})$;

$V_{j\kappa}=(\widetilde{x}_{j\kappa},\widetilde{y}_{j\kappa}),~{}\kappa=12,23,31$
are the three vertices of $T_{j}$;

$T_{jk},~{}k=1,2,3$ are the neighboring triangles that share a common side
with $T_{j}$;

$\ell_{jk}$ is the length of the common side of $T_{j}$ and $T_{jk}$, and
$M_{jk}$ is its midpoint;

$\bm{n}_{jk}:=(\cos(\theta_{jk}),\sin(\theta_{jk}))^{\top}$ is the outer unit
normal to the $k$th side of $T_{j}$.

Next, the bottom topography $B$ is replaced with its continuous piecewise
linear approximation $\widetilde{B}$ given by

$\left|\begin{array}[]{ccc}x-\widetilde{x}_{j12}&y-\widetilde{y}_{j12}&{\widetilde{B}}(x,y)-\widehat{B}_{j12}\\\\[2.15277pt]
\widetilde{x}_{j23}-\widetilde{x}_{j12}&\widetilde{y}_{j23}-\widetilde{y}_{j12}&\widehat{B}_{j23}-\widehat{B}_{j12}\\\\[2.15277pt]
\widetilde{x}_{j13}-\widetilde{x}_{j12}&\widetilde{y}_{j13}-\widetilde{y}_{j12}&\widehat{B}_{j13}-\widehat{B}_{j12}\end{array}\right|=0,\quad(x,y)\in
T_{j},$

where, in the case of continuous bottom topography,
$\widehat{B}_{j\kappa}:=B(V_{j\kappa})~{}\kappa=12,23,31$. Then, denote:

$B_{jk}:=\widetilde{B}(M_{jk}),\quad
B_{j}:=\widetilde{B}(x_{j},y_{j})=\frac{1}{3}(\widehat{B}_{j12}+\widehat{B}_{j23}+\widehat{B}_{j13}).$

At time $t$, define by $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}(t)$ the
approximation of the cell averages of the solution,

$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}(t)\approx\frac{1}{|T_{j}|}\iint\limits_{T_{j}}\bm{U}(x,y,t)\,dxdy.$

From now on, in order to shorten the formula, we suppress the time-dependence
in the algorithm. All indexed quantities used in the following formula will be
computed at time $t$. Then, as shown in [29, 4], the semi-discrete second-
order central-upwind scheme for the Saint-Venant system (1.1a)-(1.1d) on
triangular grids is written as the following system of ODEs,

$\displaystyle\frac{d\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}}{dt}=$
$\displaystyle-\frac{1}{|T_{j}|}\big{[}\bm{H}_{j1}+\bm{H}_{j2}+\bm{H}_{j3}\big{]}+\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{S}$\kern-0.50003pt}}}_{j},$ (2.3)

where the numerical fluxes through the edges of the triangular cell $T_{j}$
are

$\displaystyle\bm{H}_{jk}=$
$\displaystyle\frac{\ell_{jk}\cos(\theta_{jk})}{a_{jk}^{\rm in}+a_{jk}^{\rm
out}}\Big{[}{a_{jk}^{\rm in}\bm{F}(\bm{U}_{jk}(M_{jk}),B_{jk})+a_{jk}^{\rm
out}\bm{F}(\bm{U}_{j}(M_{jk}),B_{jk})}\Big{]}$ (2.4) $\displaystyle+$
$\displaystyle\frac{\ell_{jk}\sin(\theta_{jk})}{a_{jk}^{\rm in}+a_{jk}^{\rm
out}}\Big{[}{a_{jk}^{\rm in}\bm{G}(\bm{U}_{jk}(M_{jk}),B_{jk})+a_{jk}^{\rm
out}\bm{G}(\bm{U}_{j}(M_{jk}),B_{jk})}\Big{]}$ (2.5) $\displaystyle-$
$\displaystyle\ell_{jk}\frac{a_{jk}^{\rm in}a_{jk}^{\rm out}}{a_{jk}^{\rm
in}+a_{jk}^{\rm
out}}\big{[}\bm{U}_{jk}(M_{jk})-\bm{U}_{j}(M_{jk})\big{]},\quad k=1,2,3.$
(2.6)

Here, $\bm{U}_{j}(M_{jk})$ and $\bm{U}_{jk}(M_{jk})$ are the reconstructed
point values of $\bm{U}$ at the middle points of the edges $M_{jk}$. To obtain
these values, first, a piecewise linear reconstruction of the variables
$\bm{\Upsilon}:=(w,u,v,\rho)^{\top}$ is computed as,

$\widetilde{\bm{\Upsilon}}(x,y)=\sum\limits_{j}\bm{\Upsilon}_{j}(x,y){\Large{\mbox{$\chi$}}}_{T_{j}},\quad\bm{\Upsilon}_{j}(x,y):=\bm{\Upsilon}_{j}+(\widehat{\bm{\Upsilon}}_{x})_{j}(x-x_{j})+(\widehat{\bm{\Upsilon}}_{y})_{j}(y-y_{j}),$
(2.7)

where ${\Large{\mbox{$\chi$}}}_{T_{j}}$ is the characteristic function of the
cell $T_{j}$, $\bm{\Upsilon}_{j}$ are the point values of $\bm{\Upsilon}$ at
the cell centers, and $(\widehat{\bm{\Upsilon}}_{x})_{j}$ and
$(\widehat{\bm{\Upsilon}}_{y})_{j}$ are the limited partial derivatives (see
Section 3 in [29] for more details of the computation). In order to compute
the cell center point values of the density,
$\rho_{j}\approx\rho(x_{j},y_{j},t)$ in cell $T_{j}$ and velocities
$u_{j}\approx u(x_{j},y_{j},t)$ and $v_{j}\approx v(x_{j},y_{j},t)$, we use
the desingularization procedure presented in [29]. After that, the second and
third components of the point values $\bm{U}_{j}(M_{jk})$ and
$\bm{U}_{jk}(M_{jk})$ are obtained from, $\bm{\Upsilon}_{j}(M_{jk})$ and
$\bm{\Upsilon}_{jk}(M_{jk})$,

$\displaystyle(hu)_{j}(M_{jk})=(w_{j}(M_{jk})-B_{jk})\,u_{j}(M_{jk}),$
$\displaystyle(hu)_{jk}(M_{jk})=(w_{jk}(M_{jk})-B_{jk})\,u_{jk}(M_{jk}),$
$\displaystyle(hv)_{j}(M_{jk})=(w_{j}(M_{jk})-B_{jk})\,v_{j}(M_{jk}),$
$\displaystyle(hv)_{jk}(M_{jk})=(w_{jk}(M_{jk})-B_{jk})\,v_{jk}(M_{jk}).$
$\displaystyle(h\rho)_{j}(M_{jk})=(w_{j}(M_{jk})-B_{jk})\,\rho_{j}(M_{jk}),$
$\displaystyle(h\rho)_{jk}(M_{jk})=(w_{jk}(M_{jk})-B_{jk})\,\rho_{jk}(M_{jk}).$

However, in some numerical examples, the oscillation may appear when we use
the piecewise linear approximation (2.7) to obtain the point values in some
cells. It may occur due to the appearance of local extrema values at middle
points or as a consequence of the density discontinuity. To prevent this
oscillation, we will use some techniques that will be discussed later in
Sections 2.1 and 2.2.

In (2.6), $a_{jk}^{\rm in}$ and $a_{jk}^{\rm out}$ are the one-sided local
speeds of propagation in the directions $\pm\bm{n}_{jk}$. These speeds are
related to the largest and smallest eigenvalues of the Jacobian matrix
$J_{jk}=\cos(\theta_{jk})\,\frac{\partial{\bm{F}}}{\partial{\bm{U}}}+\sin(\theta_{jk})\,\frac{\partial{\bm{G}}}{\partial{\bm{U}}}$,
denoted by $\lambda_{+}[J_{jk}]$ and $\lambda_{-}[J_{jk}]$, respectively, and
are defined by

$\displaystyle a^{\rm
in}_{jk}=-\min\\{\lambda_{-}[J_{jk}(\bm{U}_{j}(M_{jk}))],\,\lambda_{-}[J_{jk}(\bm{U}_{jk}(M_{jk})],\,0\\},$
(2.8) $\displaystyle a^{\rm
out}_{jk}=\max\\{\lambda_{+}[J_{jk}(\bm{U}_{j}(M_{jk}))],\,\lambda_{+}[J_{jk}(\bm{U}_{jk}(M_{jk})],\,0\\},$

where

$\displaystyle\lambda_{\pm}[J_{jk}(\bm{U}_{j}(M_{jk}))]=\cos(\theta_{jk})u_{j}(M_{jk})+\sin(\theta_{jk})v_{j}(M_{jk})\pm\sqrt{\frac{g}{\rho_{0}}h_{j}(M_{jk})\rho_{j}(M_{jk})},$
$\displaystyle\lambda_{\pm}[J_{jk}(\bm{U}_{jk}(M_{jk}))]=\cos(\theta_{jk})u_{jk}(M_{jk})+\sin(\theta_{jk})v_{jk}(M_{jk})\pm\sqrt{\frac{g}{\rho_{0}}h_{jk}(M_{jk})\rho_{jk}(M_{jk})}.$

###### Remark 2.1

In order to avoid division by 0 (or by a very small positive number), the
numerical flux (2.6) is replaced with

$\displaystyle\bm{H}_{jk}=$
$\displaystyle\frac{\ell_{jk}\cos(\theta_{jk})}{2}\left[\bm{F}(\bm{U}_{jk}(M_{jk}),B_{jk})+\bm{F}(\bm{U}_{j}(M_{jk}),B_{jk})\right]$
$\displaystyle+$
$\displaystyle\frac{\ell_{jk}\sin(\theta_{jk})}{2}\left[\bm{G}(\bm{U}_{jk}(M_{jk}),B_{jk})+\bm{G}(\bm{U}_{j}(M_{jk}),B_{jk})\right]$

wherever $a_{jk}^{\rm in}+a_{jk}^{\rm out}<\sigma$. In all of the reported
numerical examples in Section 5, we have taken $\sigma=10^{-6}$.

Finally, the cell average of the source term $\bm{S}_{j}$ in (2.3),

$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{S}$\kern-0.50003pt}}}_{j}(t)\approx\frac{1}{|T_{j}|}\iint\limits_{T_{j}}\bm{S}\big{(}\bm{U}(x,y,t),B(x,y)\big{)}\,dxdy,$

has to be discretized in a well-balanced manner presented later in Section 2.4

###  2.1 The Limiter for Piecewise Linear Approximation

As discussed above, we may observe the oscillation when using piecewise linear
approximation (2.7) due to the appearance of local extrema at midpoints
$M_{jk}$, see [29]. Hence, in order to maintain the numerical stability of the
scheme, the following monotonicity condition must be satisfied.

$\min\big{(}\bm{\Upsilon}^{(i)}_{j},\bm{\Upsilon}^{(i)}_{jk}\big{)}\leq\bm{\Upsilon}^{(i)}_{j}(M_{jk})\leq\max\big{(}\bm{\Upsilon}^{(i)}_{j},\bm{\Upsilon}^{(i)}_{jk}\big{)},\quad
k=1,2,3,$ (2.9)

where $\bm{\Upsilon}^{(i)}_{j}(M_{jk})$ is the $i$-th component of the point
values at midpoint $M_{jk}$ which is obtained by the linear reconstruction
(2.7). Here, $\bm{\Upsilon}^{(i)}_{j}$ and $\bm{\Upsilon}^{(i)}_{jk}$ are the
center point value in cell $T_{j}$ and the neighboring cell $T_{jk}$, see
[19]. To ensure that the reconstructed point value
$\bm{\Upsilon}^{(i)}_{j}(M_{jk})$ is between two cell averages
$\bm{\Upsilon}^{(i)}_{j}$ and $\bm{\Upsilon}^{(i)}_{jk}$, in [29, 3], the
gradients $\widehat{\bm{\Upsilon}}_{x}$ and $\widehat{\bm{\Upsilon}}_{y}$ are
set to zero for cells violating at least one of the inequalities (2.9).
However, using zero gradients may reduce the convergence rate in numerical
simulation. Hence, in our research, instead of zero gradients, we consider the
idea of the positivity-preserving limiter developed in [32, 43] to correct the
piecewise linear reconstruction. The details for this correction are presented
as follows.

Consider an arbitrary cell $T_{j}$ with the polynomial
$\bm{\Upsilon}^{(i)}_{j}(x,y)$ that does not satisfy at least one of local
extrema conditions (2.9). In that case, we will replace
$\bm{\Upsilon}^{(i)}_{j}(x,y)$ by

$(\bm{\Upsilon}^{(i)}_{j})^{*}(x,y)=\theta(\bm{\Upsilon}^{(i)}_{j}(x,y)-\bm{\Upsilon}^{(i)}_{j})+\bm{\Upsilon}^{(i)}_{j},$
(2.10)

where $\theta\in[0,1]$ such that $(\bm{\Upsilon}^{(i)}_{j})^{*}(x,y)$
satisfies the inequalities (2.9) for all $k=1,2,3$. We then need to solve for
the constant $\theta$. Plugging (2.10) into the inequalities (2.9), we have

$\theta\in\left[\min(\alpha_{jk},\beta_{jk}),\max(\alpha_{jk},\beta_{jk})\right],\quad
k=1,2,3$ (2.11)

or

$\theta\in[\max_{k}(\min(\alpha_{jk},\beta_{jk})),\min_{k}(\max(\alpha_{jk},\beta_{jk}))]$
(2.12)

where
$\alpha_{jk}=\dfrac{\min\big{(}\bm{\Upsilon}^{(i)}_{j},\bm{\Upsilon}^{(i)}_{jk}\big{)}-\bm{\Upsilon}^{(i)}_{j}}{\bm{\Upsilon}^{(i)}_{j}(M_{jk})-\bm{\Upsilon}^{(i)}_{j}},\quad\beta_{jk}=\dfrac{\max\big{(}\bm{\Upsilon}^{(i)}_{j},\bm{\Upsilon}^{(i)}_{jk}\big{)}-\bm{\Upsilon}^{(i)}_{j}}{\bm{\Upsilon}^{(i)}_{j}(M_{jk})-\bm{\Upsilon}^{(i)}_{j}}.$

Hence we can take

$\theta=\min(\min_{k}(\max(\alpha_{jk},\beta_{jk})),1),$ (2.13)

Note that $\theta\geq 0$ since $\theta=0$, which corresponds to the constant
approximation for $(\bm{\Upsilon}^{(i)}_{j})^{*}(x,y)$, is one solution of the
inequalities (2.9). Hence, with this optimized limiter, the point values stay
within the given range. In addition, the approximated values of water depth
and density at middle points are also positive as
$\bm{\Upsilon}^{(i)}_{j}(M_{jk})\geq\min\big{(}\bm{\Upsilon}^{(i)}_{j},\bm{\Upsilon}^{(i)}_{jk}\big{)}>0$.
Therefore, the stability is ensured without using zero-gradients.

###  2.2 Riemann Solver for Mixed Cells

As mentioned in Section 1 and [9], a good scheme should not develop spurious
pressure oscillations in the neighborhood of density jumps. The oscillations
appear when we numerically solve the compressible multifluid systems by using
conventional Godunov-type finite-volume methods. We define the cells which
contains the curves of density discontinuity as “mixed” cells. Note that the
quantities in mixed cells are calculated as an artificial numerical mixture of
two different fluids. Hence, the cell averages of mixed cells may have no or
very little physical sense and become ’unreliable’ [9]. Consequently, using
the cell averages in mixed cells to obtain the linear approximation (2.7) in
the neighborhood of the interface will lead to unexpected pressure
oscillations. In [9], to eliminate the oscillation, instead of the linear
approximation (2.7), the point values in the “mixed” cells are calculated by
using the approximated solution of 1-D Riemann problems. Namely, according to
this approach, the 1-D Riemann Solver is applied in the $x$-direction to
approximate the point values on the vertical boundaries of each rectangular
mixed cell. Similarly, for the point values on horizontal sides of the
rectangular, 1-D Riemann problem is considered in $y$-direction. This leads to
the idea of applying Riemann Solver approach in the direction of the normal
vector, so called “normal direction”, to get the middle point values on each
side of the mixed cells in triangular meshes. We can understand this idea as
rotating the reference coordinate. Recently, a method using rotated coordinate
frame to solve multi-dimensional Riemann problems has been proposed in [1]. In
such approach, the 2-D Riemann Problem is converted to 1-D problem in a
particular direction in order to easily obtain the corresponding solution in
2-D. Therefore, in our research, we will calculate the point values in
triangles based on the intermediate states of 1-D Riemann problem in the
normal direction as follows.

Assume that triangle $T_{j}$ is a mixed cell which has the outward unit normal
vector $\mathbf{n}_{jk}=(\cos(\theta_{jk}),\sin(\theta_{jk}))$ on side $k$. We
will consider a new reference frame $(x^{\prime},y^{\prime})$ such that the
new horizontal axis is in the direction of the outward normal vector
$\mathbf{n}_{jk}$ and the origin $(0,0)$ is at the middle point $M_{jk}$ of
side $k$ as shown in Fig. 2.2.

Figure 2.2: An example of the rotated coordinates $(x^{\prime},y^{\prime})$
where $x^{\prime}$-axis is in the direction of normal vector $\mathbf{n}_{jk}$
and the new origin is at midpoint $M_{jk}$ of side $k$.

Suppose $T_{L}$ and $T_{R}$ to be the two closest single-fluid cells to
$T_{j}$ such that $T_{L}$ and $T_{R}$ stay on different sides of the $k$-th
edge of cell $T_{j}$. We assume that $T_{R}$ stays on the side where the
outward normal vector $\mathbf{n}_{jk}$ is pointing. Let
$\bm{\Upsilon}_{L}=(w_{L},u_{L},v_{L},\rho_{L})$ and
$\bm{\Upsilon}_{R}=(w_{R},u_{R},v_{R},\rho_{R})$ be the center point values of
single-fluid cells $T_{L}$ and $T_{R}$. Note that $(u_{L},v_{L})$ and
$(u_{R},v_{R})$ are the values of velocities in the Cartesian coordinate
system $(x,y)$. The projections of velocities $(u_{L},v_{L})$ and
$(u_{R},v_{R})$ onto the rotated frame $(x^{\prime},y^{\prime})$ are

$u^{\prime}_{L}=u_{L}\cos(\theta_{jk})+v_{L}\sin(\theta_{jk}),\quad
v^{\prime}_{L}=-u_{L}\sin(\theta_{jk})+v_{L}\cos(\theta_{jk}),$
$u^{\prime}_{R}=u_{R}\cos(\theta_{jk})+v_{R}\sin(\theta_{jk}),\quad
v^{\prime}_{R}=-u_{R}\sin(\theta_{jk})+v_{R}\cos(\theta_{jk}).$

Similar to [1, 9], to compute the point value at midpoint $M_{jk}$ on side $k$
of triangle $T_{j}$, we have to solve the following 1-D Riemann problems
between states $\bm{\Upsilon}_{L}$ an $\bm{\Upsilon}_{R}$ in direction
$\mathbf{n}_{jk}$.

$\begin{cases}w_{t}+(hu)_{x^{\prime}}=0\\\
(hu)_{t}+\Big{(}\dfrac{(hu)^{2}}{w-B}+\frac{g}{2}\frac{\rho}{\rho_{0}}(w-B)\Big{)}_{x^{\prime}}=-g\dfrac{\rho}{\rho_{0}}(w-B)B_{x^{\prime}},\\\
(hv)_{t}+\left(\dfrac{(hu)(hv)}{w-B}\right)_{x^{\prime}}=0,\\\
(h\frac{\rho}{\rho_{0}})_{t}+\left(uh\dfrac{\rho}{\rho_{0}}\right)_{x^{\prime}}=0,\end{cases}$
(2.14)

subject to the following initial condition

$(w,u,v,\rho,B)(x^{\prime},0)=\begin{cases}(\bm{\Upsilon}_{L})^{\prime}:=(w_{L},u^{\prime}_{L},v^{\prime}_{L},\rho_{L},B_{L})\quad\mbox{if}\quad
x^{\prime}<0,\\\
(\bm{\Upsilon}_{R})^{\prime}:=(w_{R},u^{\prime}_{R},v^{\prime}_{R},\rho_{R},B_{R})\quad\mbox{if}\quad
x^{\prime}>0.\end{cases}$ (2.15)

The details of the approximate Riemann solver for the above Riemann problem
can be seen in [9]. Suppose
$(\bm{\Upsilon}^{*}_{L})^{\prime}:=(w^{*}_{L},(u^{*}_{L})^{\prime},(v^{*}_{L})^{\prime},\rho^{*}_{L})$
and
$(\bm{\Upsilon}^{*}_{R})^{\prime}:=(w^{*}_{R},(u^{*}_{R})^{\prime},(v^{*}_{R})^{\prime},\rho^{*}_{R})$
are the intermediate states calculated by Riemann solver for (2.14)-(2.15) in
$(x^{\prime},y^{\prime})$ coordinates. We recall that in [9], the point values
at midpoints on left and right boundaries of a rectangular mixed cell are
obtained respectively based on the left and right Riemann intermediate
states,$(\bm{\Upsilon}^{*}_{L})^{\prime}$ and
$(\bm{\Upsilon}^{*}_{R})^{\prime}$. Here, in triangular grids, we can choose
either $(\bm{\Upsilon}^{*}_{L})^{\prime}$ or
$(\bm{\Upsilon}^{*}_{R})^{\prime}$ to calculate the point values at the middle
point $M_{jk}$ of cell $T_{j}$. Note that $T_{L}$ and $T_{R}$ are both single-
fluid cells nearby mixed cell $T_{j}$. However, if the neighboring cell
$T_{jk}$ shown in Fig. 2.2 is a single-fluid cell, then $T_{R}\equiv T_{jk}$
and $M_{jk}\in T_{R}$, which means $T_{R}$ is closer to $M_{jk}$ compared to
$T_{L}$. Hence, in our work, we select the right intermediate values
$(\bm{\Upsilon}^{*}_{R})^{\prime}$ and the information from $T_{R}$ to compute
the point values $\bm{\Upsilon}_{j}(M_{jk})$. From [9], we first calculate the
sound of speed $c^{*}_{R}=h^{2}\dfrac{g\rho}{2\rho_{0}}$. If
$h^{*}_{R}>0,\rho^{*}_{R}>0$ and $(u^{*})^{\prime}_{R}-c^{*}_{R}<0$, then we
replace the piecwise linear approximation in (2.7) by
$\bm{\Upsilon}_{j}(M_{jk})=(\bm{\Upsilon}^{*}_{R})^{\prime}$, otherwise
$\bm{\Upsilon}_{j}(M_{jk})=(\bm{\Upsilon}_{R})^{\prime}$ (more details of the
point values correction can be seen in [9]). Finally, we convert the computed
solution back to original Cartesian coordinate $(x,y)$ as

$\bm{\Upsilon}_{j}(M_{jk})=\begin{cases}\bm{\Upsilon}^{*}_{R},\quad&\mbox{if}\quad
h^{*}_{R}>0,\rho^{*}_{R}>0,(u^{*})^{\prime}_{R}-c^{*}_{R}<0,\\\
\bm{\Upsilon}_{R},\quad&\mbox{otherwise},\end{cases}$ (2.16)

where
$\bm{\Upsilon}^{*}_{R}:=(w^{*}_{R},(u^{*}_{R})^{\prime}\cos(\theta_{jk})-(v^{*}_{R})^{\prime}\sin(\theta_{jk}),(u^{*}_{R})^{\prime}\sin(\theta_{jk})+(v^{*}_{R})^{\prime}\cos(\theta_{jk}),\rho^{*}_{R})$.

From [9], in “lake at rest” solutions (1.2) and (1.3), we have
$\bm{\Upsilon}^{*}_{R}=\bm{\Upsilon}_{R}$ and the steady state solutions are
then preserved by using the discretization of the source term presented in
Section 2.4.

###  2.3 Positivity-preserving condition

Next, we will discuss the restriction of time step to preserve the positivity
of water depth $h$ and the variable $h\rho$. Consider the Forward Euler (FE)
equation

$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt${U}$\kern-0.50003pt}}}_{j}^{n+1}=\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt${U}$\kern-0.50003pt}}}_{j}^{n}-\frac{1}{|T_{j}|}\sum_{k=1}^{3}\Delta
t_{jk}{H}_{jk}+\Delta t\,\bar{{S}}_{j}.$

From [3], the time step condition to guarantee the positivity of water height
$h$ is

$\Delta t<\dfrac{1}{6a}\min_{jk}(r_{jk}).$ (2.17)

One can clearly see that the time step size (2.17) will also achieve
nonnegative water height $h$ in the density shallow water model (1.1a)-(1.1d),
see [3, 29]. We now use the idea from [3] to find the CFL-type condition for
the positivity-preserving property of the last variable $h\rho$. To this end,
we apply the forward Euler discretization to the last component of the scheme.

$\displaystyle(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}})^{n+1}_{j}=(\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}})^{n}_{j}-$
$\displaystyle\dfrac{\Delta
t}{|T_{j}|}\sum_{k=1}^{3}\frac{\ell_{jk}\cos(\rho_{jk})}{a_{jk}^{\rm
in}+a_{jk}^{\rm out}}\Big{[}a_{jk}^{\rm in}(h\rho u)_{jk}(M_{jk})+a_{jk}^{\rm
out}(h\rho u)_{j}(M_{jk})\Big{]}$ (2.18) $\displaystyle-$
$\displaystyle\dfrac{\Delta
t}{|T_{j}|}\sum_{k=1}^{3}\frac{\ell_{jk}\sin(\rho_{jk})}{a_{jk}^{\rm
in}+a_{jk}^{\rm out}}\Big{[}a_{jk}^{\rm in}(h\rho v)_{jk}(M_{jk})+a_{jk}^{\rm
out}(h\rho v)_{j}(M_{jk})\Big{]}$ $\displaystyle+$ $\displaystyle\dfrac{\Delta
t}{|T_{j}|}\sum_{k=1}^{3}\ell_{jk}\frac{a_{jk}^{\rm in}a_{jk}^{\rm
out}}{a_{jk}^{\rm in}+a_{jk}^{\rm
out}}\big{[}(h\rho)_{jk}(M_{jk})-(h\rho)_{j}(M_{jk})\big{]}.$

Note that

$(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}})^{n}_{j}=h_{j}\rho_{j}=\left(\dfrac{1}{3}\sum_{m=1}^{3}h_{j}(M_{jm})\right)\left(\dfrac{1}{3}\sum_{s=1}^{3}\rho_{j}(M_{js})\right)=\dfrac{1}{9}\sum_{m,s}h_{j}(M_{jm})\rho_{j}(M_{js}).$
(2.19)

Plug (2.19) to (2.18), we have

$\displaystyle(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}})^{n+1}_{j}=$
$\displaystyle\dfrac{1}{9}\sum_{m,s}h_{j}(M_{jm})\rho_{j}(M_{js})$
$\displaystyle-\dfrac{\Delta
t}{|T_{j}|}\sum_{k=1}^{3}\frac{\ell_{jk}\cos(\rho_{jk})}{a_{jk}^{\rm
in}+a_{jk}^{\rm out}}\Big{[}a_{jk}^{\rm
in}h_{jk}(M_{jk})\rho_{jk}(M_{jk})u_{jk}(M_{jk})+a_{jk}^{\rm
out}h_{j}(M_{jk})\rho_{j}(M_{jk})u_{j}(M_{jk})\Big{]}$
$\displaystyle-\dfrac{\Delta
t}{|T_{j}|}\sum_{k=1}^{3}\frac{\ell_{jk}\sin(\rho_{jk})}{a_{jk}^{\rm
in}+a_{jk}^{\rm out}}\Big{[}a_{jk}^{\rm
in}h_{jk}(M_{jk})\rho_{jk}(M_{jk})v_{jk}(M_{jk})+a_{jk}^{\rm
out}h_{j}(M_{jk})\rho_{j}(M_{jk})v_{j}(M_{jk})\Big{]}$
$\displaystyle+\dfrac{\Delta
t}{|T_{j}|}\sum_{k=1}^{3}\ell_{jk}\frac{a_{jk}^{\rm in}a_{jk}^{\rm
out}}{a_{jk}^{\rm in}+a_{jk}^{\rm
out}}\big{[}h_{jk}(M_{jk})\rho_{jk}(M_{jk})-h_{j}(M_{jk})\rho_{j}(M_{jk})\big{]}$
$\displaystyle=$ $\displaystyle\dfrac{1}{9}\sum_{m\neq
s}h_{j}(M_{jm})\rho_{j}(M_{js})$
$\displaystyle+\sum_{k=1}^{3}h_{j}(M_{jk})\rho_{j}(M_{jk})\Big{[}\dfrac{1}{9}-\dfrac{\Delta
t}{|T_{j}|}\frac{\ell_{jk}a_{jk}^{\rm out}}{a_{jk}^{\rm in}+a_{jk}^{\rm
out}}\left(a_{jk}^{\rm in}+u^{\perp}_{j}(M_{jk})\right)\Big{]}$
$\displaystyle+\sum_{k=1}^{3}h_{jk}(M_{jk})\rho_{jk}(M_{jk})\dfrac{\Delta
t}{|T_{j}|}\frac{\ell_{jk}a_{jk}^{\rm in}}{a_{jk}^{\rm in}+a_{jk}^{\rm
out}}\left(a_{jk}^{\rm out}-u^{\perp}_{jk}(M_{jk})\right).$

From the definitions of the local speeds (2.8) we obtain that
$u^{\perp}_{j}(M_{jk})=\cos(\theta_{jk})u_{j}(M_{jk})+\sin(\theta_{jk})v_{j}(M_{jk})\leq
a_{jk}^{\rm out}$, $0\leq a_{jk}^{\rm in}$, and $0\leq a_{jk}^{\rm out}$.
Besides, the corrected reconstruction in Section 2.1 guarantees that
$h_{jk}(M_{jk})\geq 0$ and $\rho_{jk}(M_{jk})\geq 0$ for all $j$ and
$k=1,2,3$. Therefore, all terms in the first and third sum on the RHS of (2.3)
are nonnegative. To enforce the second sum on the RHS of (2.3) to be
nonegative, we then need:

$\dfrac{\Delta t}{|T_{j}|}\frac{\ell_{jk}a_{jk}^{\rm out}}{a_{jk}^{\rm
in}+a_{jk}^{\rm out}}\left(a_{jk}^{\rm
in}+u^{\perp}_{j}(M_{jk})\right)\leq\dfrac{\Delta
t}{|T_{j}|}\frac{\ell_{jk}a_{jk}^{\rm out}}{a_{jk}^{\rm in}+a_{jk}^{\rm
out}}\left(a_{jk}^{\rm in}+a_{jk}^{\rm out}\right)=\dfrac{\Delta
t}{|T_{j}|}l_{jk}a_{jk}^{\rm out}<\dfrac{1}{9},\quad\forall j\mbox{ and
}k=1,2,3.$

Hence we conclude that all terms in the second sum on the RHS of (2.3) are
also nonnegative under the following CFL-type condition.

$\Delta
t<\dfrac{1}{18a}\min_{jk}(r_{jk})=\dfrac{\min_{jk}\left(\dfrac{2|T_{j}|}{l_{jk}}\right)}{18\max_{jk}(a_{jk}^{\rm
in},a_{jk}^{\rm out})}<\dfrac{1}{9}\min_{jk}\left(\dfrac{|T_{j}|}{l_{jk}a^{\rm
out}_{jk}}\right),$ (2.20)

where $a=\max_{jk}(a_{jk}^{\rm in},a_{jk}^{\rm out})$ and
$r_{jk}=\dfrac{2|T_{j}|}{l_{jk}}$ is the $k$-th altitude of triangle $T_{j}$.
This completes the proof. This proof is still valid if one uses a higher-order
SSP ODE solver (either the Runge-Kutta or the multistep one), because such
solvers can be written as a convex combination of several forward Euler steps.

###  2.4 Well-balanced discretization of source term.

In this section, we first develop a well-balanced discretization of the source
terms, which guarantees the first type of “lake at rest” steady-state
solutions,

$w=\max\big{\\{}C,B(x,y)\big{\\}},\quad C=\rm{Const},\quad\rho=P\equiv
Const,\quad u\equiv v\equiv 0,$ (2.21)

are exactly preserved by the resulting central-upwind scheme. This means that
the source discretization $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\mathbf{S}$\kern-0.50003pt}}}_{j}$ should
exactly balance the numerical fluxes so that the right-hand side (RHS) of
(1.1a)-(1.1d) vanishes at “lake at rest” steady states.

To this end, we substitute the “lake at rest” state (2.21) into FE scheme of
the system (it is still correct with the adaptive Runge-Kutta scheme) and
conclude that a well-balanced quadrature for $\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\mathbf{S}$\kern-0.50003pt}}}_{j}$ should
satisfy the following two conditions:

$-\frac{g}{2|T_{j}|}\sum_{k=1}^{3}\ell_{jk}\cos(\theta_{jk})\cdot\frac{\Delta
t^{(\rho)}_{jk}}{\Delta
t}\cdot\dfrac{P}{\rho_{0}}\big{[}C-B(M_{jk})\big{]}^{2}+\bar{S}_{j}^{\,(2)}=0,$
(2.22)

$-\frac{g}{2|T_{j}|}\sum_{k=1}^{3}\ell_{jk}\sin(\theta_{jk})\cdot\frac{\Delta
t^{(\rho)}_{jk}}{\Delta
t}\cdot\dfrac{P}{\rho_{0}}\big{[}C-B(M_{jk})\big{]}^{2}+\bar{S}_{j}^{\,(3)}=0,$
(2.23)

where

$\bar{S}_{j}^{\,(2)}\approx-\frac{g}{|T_{j}|\rho_{0}}\iint\limits_{T_{j}}P(C-B(x,y))B_{x}(x,y)\,dxdy,\quad\bar{S}_{j}^{\,(3)}\approx-\frac{g}{|T_{j}|\rho_{0}}\iint\limits_{T_{j}}P(C-B(x,y))B_{y}(x,y)\,dxdy.$

In order to derive a well-balanced quadrature, similar to [3, 29], we first
apply Green’s formula, $\iint_{T_{j}}{\rm div}\,\mathbf{{\cal
G}}\,dxdy=\int_{\partial T_{j}}\mathbf{{\cal G}}\cdot\mathbf{n}\,ds$, to the
vector field $\mathbf{{\cal G}}=(\dfrac{1}{2}\rho(x,y)(w(x,y)-B(x,y))^{2},0)$
and obtain

$\displaystyle-\iint\limits_{T_{j}}\rho(x,y)(w(x,y)-B(x,y))B_{x}(x,y)\,dxdy=$
$\displaystyle\sum_{k=1}^{3}\int\limits_{(\partial{T_{j}})_{k}}\rho(x,y)\frac{(w(x,y)-B(x,y))^{2}}{2}\cos(\theta_{jk})\,ds$
(2.24)
$\displaystyle-\iint\limits_{T_{j}}\rho(x,y)(w(x,y)-B(x,y))w_{x}(x,y)\,dxdy$
$\displaystyle-\iint\limits_{T_{j}}\rho_{x}(x,y)\frac{(w(x,y)-B(x,y))^{2}}{2}\,dxdy,$

where $(\partial{T_{j}})_{k}$ is the $k$-th side of the triangle $T_{j}$,
$k=1,2,3$. The double integrals are approximated using the trapezoidal rule.
This results in the following quadrature for $\bar{S}^{\,(2)}_{j}$:

$\displaystyle\bar{S}^{\,(2)}_{j}$
$\displaystyle=\frac{g}{2|T_{j}|}\sum_{k=1}^{3}\ell_{jk}\cos(\theta_{jk})\cdot\frac{\Delta
t^{(\rho)}_{jk}}{\Delta
t}\cdot\dfrac{\rho(M_{jk})}{\rho_{0}}\big{[}w(M_{jk})-B(M_{jk})\big{]}^{2}$
(2.25)
$\displaystyle-\frac{g}{3\rho_{0}}\Big{[}\rho_{j12}(w_{j12}-\widehat{B}_{j12})\,w_{x}(V_{j12})+\rho_{j23}(w_{j23}-\widehat{B}_{j23})\,w_{x}(V_{j23})+\rho_{j13}(w_{j13}-\widehat{B}_{j13})\,w_{x}(V_{j13})\Big{]}.$
$\displaystyle-\frac{g}{6\rho_{0}}\Big{[}\rho_{x}(V_{j12})(w_{j12}-\widehat{B}_{j12})^{2}\,+\rho_{x}(V_{j23})(w_{j23}-\widehat{B}_{j23})^{2}\,+\rho_{x}(V_{j13})(w_{j13}-\widehat{B}_{j13})^{2}\,\Big{]}.$

A similar quadrature for $\bar{S}^{\,(3)}_{j}$ is

$\displaystyle\bar{S}^{\,(3)}_{j}$
$\displaystyle=\frac{g}{2|T_{j}|}\sum_{k=1}^{3}\ell_{jk}\sin(\theta_{jk})\cdot\frac{\Delta
t^{(\rho)}_{jk}}{\Delta
t}\cdot\dfrac{\rho(M_{jk})}{\rho_{0}}\big{[}w(M_{jk})-B(M_{jk})\big{]}^{2}$
(2.26)
$\displaystyle-\frac{g}{3\rho_{0}}\Big{[}\rho_{j12}(w_{j12}-\widehat{B}_{j12})\,w_{y}(V_{j12})+\rho_{j23}(w_{j23}-\widehat{B}_{j23})\,w_{y}(V_{j23})+\rho_{j13}(w_{j13}-\widehat{B}_{j13})\,w_{y}(V_{j13})\Big{]}.$
$\displaystyle-\frac{g}{6\rho_{0}}\Big{[}\rho_{y}(V_{j12})(w_{j12}-\widehat{B}_{j12})^{2}\,+\rho_{y}(V_{j23})(w_{j23}-\widehat{B}_{j23})^{2}\,+\rho_{y}(V_{j13})(w_{j13}-\widehat{B}_{j13})^{2}\,\Big{]}.$

Notice that the piecewise linear reconstruction procedure ensures that at the
steady state (2.21), $\nabla\rho(V_{j\kappa})=0,u_{x}=v_{y}=0$ and $\nabla
w(V_{j\kappa})=0$ throughout the entire computational domain. This implies
that
$(w_{j\kappa}-\widehat{B}_{j\kappa})\,w_{x}(V_{j\kappa})\equiv(w_{j\kappa}-\widehat{B}_{j\kappa})\,w_{y}(V_{j\kappa})\equiv
0$ and
$\rho_{x}(V_{j\kappa})(w_{j\kappa}-\widehat{B}_{j\kappa})^{2}\equiv\rho_{y}(V_{j\kappa})(w_{j\kappa}-\widehat{B}_{j\kappa})^{2}\equiv
0$. Therefore, the quadratures (2.25) and (2.26) satisfy the desired well-
balanced requirements (2.22) and (2.23).

Next, we consider the “lake at rest” situation

$B\equiv\rm{Const},\quad h^{2}\rho=Q\equiv Const,\quad u\equiv v\equiv 0.$
(2.27)

Note that in the region occupied by only one fluid, the density $\rho$ is a
constant and water surface $w$ is then also a constant based on (2.27). Hence,
in single-fluid cells, the steady state (2.27) is equivalent to the solution
(2.21) which is maintained by the discretization of source term (2.25)-(2.26).
In addition, in mixed cells, the point values are obtained by the data from
nearby single-fluid cells thanks to the Riemann Solver approximation presented
in Section 2.2, see also [9]. Therefore, the discretization of source term
(2.25)-(2.26), is also capable of preserving the steady state solution (2.27).

##  3 Interface Tracking

To achieve an efficient adaptive scheme for multifluid flows, it is essential
to accurately capture the curve where two types of fluid join simultaneously
with the flow field evolution. The interface tracking is important in
preventing excessive numerical diffusion of variable density, see Section 3.3.
We also need the location of density jumps to exactly update the cell averages
in the adaptive mesh reconstruction, Section 4. Therefore, we desire a simple
and effective interface reconstruction for the system (1.1a)-(1.1d). Over the
last two decades, many methods have been proposed for this purpose such as the
level set method [31, 38], the volume of fluid method [17, 34], and the front
tracking method [42, 41]. Among various versions of interface tracking, we
have considered the approach described in [45, 44] for our work due to its
simplicity, robustness, and high efficiency. In particular, the interface is
obtained by using both the level set and the volume fraction functions. The
details of the interface reconstruction is described as follows.

###  3.1 Mixed Cell Detecting by Level Set Function

In the first part of this section, we will discuss the method of detecting the
mixed cells in the triangular grids. For this end, we consider the level set
function $\phi$ which is defined such that it is positive in one fluid, is
negative in the other fluid, and has zero value at the interface. Very often,
the level set value of each grid point is initialized by the signed distance
from that point to the curve of density discontinuities, see [29, 31, 38, 45].
As discussed in [29, 31, 38, 45], level set function is evolved by the
velocity $(u,v)$ of the flow field as

$\phi_{t}+u\phi_{x}+v\phi_{y}=0.$ (3.1)

The equation (3.1) can be rewritten in conservation form as follows.

$\phi_{t}+(u\phi)_{x}+(v\phi)_{y}=(u_{x}+v_{y})\phi.$ (3.2)

We have used the central-upwind scheme presented in Section 2 to solve for
$\phi$, where the point values $\phi(M_{jk}),k=1,2,3$ in triangle $T_{j}$ are
approximated by the piecewise linear reconstruction (2.7). The integral of the
source term in (3.2) can be approximated by midpoint rule. Note that using the
central-upwind scheme to solve (3.1) only give us the cell averages of $\phi$
in each cell and does not point out which cells contain the interface
$\phi=0$. We have to provide a numerical method for detecting mixed cells.
However, a triangle is a single-fluid cell, also called as “reliable” cell, if
the level set values at its vertices are either all positive or all negative.
Hence, we will locate each node in the grid based on its point value of level
set as presented below.

Without the loss of generality, we assume that $\phi(x,y)$ is positive in
fluid 1 and $\phi(x,y)$ is negative in fluid 2. The point value of level set
function $\phi^{*}_{j\kappa}$ at each vertex $V_{j\kappa}$ of cell $T_{j}$ is
approximated by extrapolation [9]

$\phi^{*}_{j\kappa}=\dfrac{\displaystyle\sum_{i=0}^{m^{*}_{\kappa}}c_{i}\phi^{i}_{j\kappa}}{\displaystyle\sum_{i=0}^{m^{*}_{\kappa}}c_{i}},$
(3.3)

where $\phi^{i}_{j\kappa},$ is the cell averages of level set function in
cells $T^{i}_{j\kappa}$ which have common vertex at $V_{j\kappa}$, and the
weight $c_{i}$ is inversely proportional to the distance between the center of
cells $T^{i}_{j\kappa},i=0,1,...,m^{*}_{\kappa}$ and vertex $V_{j\kappa}$. If
$\phi^{*}_{j\kappa}>0$, we have the vertex $V_{j\kappa}$ staying in fluid 1.
Otherwise, if $\phi^{*}_{j\kappa}<0$, the vertex $V_{j\kappa}$ is in fluid 2.
Finally, based on the physical meaning of the level set function, a triangle
$T_{j}$ is naturally marked as a mixed cell if it has two vertices staying in
different fluids.

###  3.2 Interface Reconstruction

In most works of interface tracking, see [31, 38, 17, 34, 42, 41], the curve
of density discontinuity in a mixed cell is approximated as a linear line
segment of the form $\mathbf{n}\cdot x=\alpha$, where $\mathbf{n}$ is the
normal vector of the interface, $x$ is the location of a point, and $\alpha$
is the line constant. A variety of methods for computing normal vector
$\mathbf{n}$ and parameter $\alpha$ have been briefly introduced and compared
in [45, 33]. Based on the advantages and disadvantages of conventional
interface tracking methods, an innovative approach has been proposed in [45].
The approach employs both the level set and volume of fluid methods to denote
the density segment in each flagged cell by its two endpoints. Namely, the
interface normal vector is calculated from the level set function while the
exact location of two endpoints of the segment is determined by enforcing mass
conservation based on the volume fraction. In our work, we will apply this
interface tracking method in the adaptive triangular mesh due to its
simplicity, accuracy, and versatility. The details of this method can be seen
in [45]. In the following, we will briefly present the main steps of the
interface reconstruction.

Once the mixed cells are detected, we will first compute the interface normal
vector in each flagged cell based on the level set function and a least square
problem. Namely, for each mixed cell $T_{j}$ with center at $(x_{j},y_{j})$,
we consider a stencil that consists of all centers $(x_{i},y_{i}),i=1,2,...,N$
of $N$ cells that share at least a vertex with $T_{j}$ (to simplify the
computation, we do not place the vertices of cell $T_{j}$ in the stencil as in
[45]). A quadratic function for function $\phi$ of a point $(x,y)$ is then
given in the generic form

$\phi=ax^{\prime 2}+bx^{\prime}y^{\prime}+cy^{\prime
2}+dx^{\prime}+ey^{\prime}+g,$

where $(x^{\prime}:=x-x_{j},y^{\prime}:=x-x_{j})$ is the location of point
$(x,y)$ in the local coordinate frame $x^{\prime}$ and $y^{\prime}$ obtained
by shifting the original coordinate such that the new origin is at the center
of $T_{j}$. The coefficients $a,b,c,d,e,$ and $g$ are determined by using the
least squares method with the linear system

$Qs=r,$

where

$Q=\begin{bmatrix}x^{\prime 2}_{1}&x^{\prime}_{1}y^{\prime}_{1}&y^{\prime
2}_{1}&x^{\prime}_{1}&y^{\prime}_{1}&1\\\ x^{\prime
2}_{2}&x^{\prime}_{2}y^{\prime}_{2}&y^{\prime
2}_{2}&x^{\prime}_{2}&y^{\prime}_{2}&1\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\ x^{\prime
2}_{N}&x^{\prime}_{N}y^{\prime}_{N}&y^{\prime
2}_{N}&x^{\prime}_{N}&y^{\prime}_{N}&1\\\ \end{bmatrix},\quad
s=\begin{bmatrix}a\\\ b\\\ c\\\ d\\\ e\\\
g\end{bmatrix},\quad\phi=\begin{bmatrix}\phi_{1}\\\ \phi_{2}\\\ \vdots\\\
\phi_{n}\end{bmatrix},$ (3.4)

$(x^{\prime}_{i}=x_{i}-x_{j},y^{\prime}_{i}=y_{i}-y_{j}),i=1,...,N,$ are the
local coordinates of the $i-$th node in the stencil and $\phi_{i},i=1,...,N,$
is the cell average of level set function in cell $T_{i}$ corresponding to
node $i$-th. Once the coefficients are known, we compute the unit normal
vector by
$\mathbf{n}=\left(\dfrac{d}{\sqrt{d^{2}+e^{2}}},\dfrac{e}{\sqrt{d^{2}+e^{2}}}\right).$
The system (3.4) is solvable as explained in [45]. In addition, since the
level set function is continuous, the normal calculation is second-order
accurate.

Next, we continue to determine the exact coordinates of two endpoints of the
interface by using the volume of fluid approach (VOF). The VOF method has been
widely used and developed for interface tracking in [17, 34, 44, 45]. In the
VOF method, the volume fraction function, denoted by $f$, is defined as the
ratio of the volume of one fluid, called fluid 1, in each cell to the total
volume of the cell. Hence, $f$ is unity if the cell is a single-fluid cell in
fluid 1, and is zero if the cell only contains fluid 2. For mixed cells, we
have $0<f<1$. The function $f$ is advected by

$f_{t}+uf_{x}+vf_{y}=0.$ (3.5)

The main idea of the VOF method proposed in [44, 45] is to reconstruct the
interface segment by determining the coordinates of its endpoints. In
particular, the two endpoints must form a segment which has the normal vector
$\mathbf{n}$ as computed above and the interface truncates the cell with the
given volume fraction. Fig. 3.1 is an illustration for two cases of a mixed
cell $T_{j}$ where the interface $I_{1}I_{2}$ splits $T_{j}$ into two parts,
$T^{(1)}_{j}$ and $T^{(2)}_{j}$, respectively occupied by fluid 1 and fluid 2.

Figure 3.1: An example of mixed cell $T_{j}=V_{j12}V_{j23}V_{j13}$ with the
interface segment $I_{1}I_{2}$ and the normal vector $\vec{\mathbf{n}}$.

We then have

$I_{1}I_{2}\cdot\mathbf{n}=0,\quad\dfrac{T^{(1)}_{j}}{T_{j}}=\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$f$\kern-0.50003pt}}}_{j},$

where $\mathbf{n}$ is the normal vector of the interface computed by (3.4) and
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$f$\kern-0.50003pt}}}_{j}$ is the cell average
of volume fraction in mixed cell $T_{j}$ obtained by using the central-upwind
method, see Section 2, to solve (3.5). This interface reconstruction is
explicit, accurate, and capable of conserving the volume of the fluid. The
details of endpoint calculation can be seen in [45]. Finally, the
reconstructed interface will be used in the adaptive mesh reconstruction as
discussed in the following sections.

###  3.3 The Cell Averages Correction

The idea of correcting the solution in the neighborhood of the interface has
been used in numerical methods for multifluid flows [9, 45] to improve the
computed solution. This technique prevents the diffusion of density emerging
when we use the central-upwind method to solve compressible systems. In our
work, we will also perform the correction procedure such that the local
conservation is ensured based on the location of density jumps. Namely,
suppose we have two types of fluid, fluid 1 and fluid 2, in the flow. After
determining the single-fluid cells and mixed cells by using the point values
of the level set function, see Section 3.1, the cell average correction will
proceed as follows.

  * •

If a cell $T_{j}$ is a single fluid cell in fluid $i=1,2$, we correct the cell
averages of $h\rho$ by $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}^{new}_{j}=h_{j}\rho^{(i)}_{j}$,
where $\rho^{(i)}_{j}$ is the value of density in fluid $i=1,2$.
$h_{j}=\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$w$\kern-0.50003pt}}}_{j}-B_{j}$ is the cell
average of water depth in $T_{j}$. To preserve the mass of $h\rho$ in the
domain, we equally split the change of  $h\rho$ in $T_{j}$ into the nearby
mixed cells. Namely, if the neighboring cell $T_{jk},k=1,2,3$ of $T_{j}$ is a
mixed cell, we obtain the new cell average of $h\rho$ in $T_{jk}$ by

$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}^{new}_{jk}=\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}_{jk}+\dfrac{\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}_{j}-\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}^{new}_{j}}{n_{mix}},$
(3.6)

where $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}_{j}$ and
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}_{jk}$ are the old cell
averages computed by using the central-upwind scheme in triangles $T_{j}$ and
$T_{jk}$, and $n_{mix}$ is the number of mixed cells surrounding the cell
$T_{j}$.

  * •

If $T_{j}$ is a mixed cell, except the update (3.6) from its nearby single-
fluid cells, no further correction is needed.

Due to the fact that the interface normally moves from one cell to its
adjacent cells, only cells surrounding the interface are updated and the loss
of mass from this correction is negligible. In addition, for cells in the
“lake at rest” area, this correction will not change the existing solution
therefore maintaining the well-balanced properties.

##  4 Adaptive Central-Upwind Scheme

The traditional numerical methods for system (1.1a)–(1.1d) consider very fine
fixed meshes to reconstruct delicate features of the solution. However, this
can lead to high computational cost, as well as to a poor accuracy of small
scale characteristics of the problem. Therefore, we will use adaptive meshes
to improve the accuracy of the approximation at a much lower cost. In this
section, we will review an efficient and accurate adaptive central-upwind
algorithm from [13] which we adapt to the system (1.1a)-(1.1d).

###  4.1 Adaptive Central-Upwind Algorithm

From [13], the adaptive central-upwind algorithm for the system (1.1a)- (1.1d)
is described briefly by the following steps.

Step 0. At time $t=t^{0}$, generate the initial uniform grid
$\mathcal{T}^{0,0}$.

Step 1. On mesh $\mathcal{T}^{n,\mathcal{M}_{n}}$, evolve the cell averages
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}^{n}$ of the solution
from time $t^{n}$ to $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}^{n+1}$ at the next
time level $t^{n+1}$ using adaptive central-upwind scheme (4.5a)-(4.5b), see
Section 4:

  * •

At time $t^{n}$, determine the level $l=0,1,...,L$ of each triangle
$T^{n,\mathcal{M}_{n}}_{j}\in\mathcal{T}^{n,\mathcal{M}_{n}}$, (4.1), Section
4.3.

  * •

At each time level $t^{n,p}_{l},p=0,1,...,\mathcal{P}_{l}-1$, perform the
piecewise polynomial reconstruction (2.7) for single-fluid cells and apply the
Riemann Solver (2.14) on the mixed cells to calculate the point values,
Section 2, Section 3.3.

  * •

At each time level $t^{n,p}_{l},p=0,1,...,\mathcal{P}_{l}-1$, compute the one-
sided local speeds of propagation using (2.8), Section 2.

  * •

At time $t^{n}$, calculate the reference time step $\Delta t$ using (4.2),
Section 4.3.

  * •

At each time level $t^{n,p}_{l},p=0,1,...,\mathcal{P}_{l}-1$, compute the
local time step for each cell level, (4.4), Section 4.3.

  * •

At each time level $t^{n,p}_{l},p=0,1,...,\mathcal{P}_{l}-1$, compute
numerical fluxes and source term in the adaptive central-upwind scheme
(4.5a)-(4.5b), (2.6), (2.25)-(2.26), Section 2, Section 4.3.

Step 2. On mesh $\mathcal{T}^{n,\mathcal{M}_{n}}$, compute WLR error using
(4.9) in Section 4.4 and determine the refinement/de-refinement status for
each cell, Section 4.4.

Step 3. Generate the new adaptive mesh $\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$
at $t^{n+1}$, Section 4.2. This step includes coarsening of some cells,
refinement of some cells, and the appropriate projection of the cell averages
from the mesh $\mathcal{T}^{n,\mathcal{M}_{n}}$ at $t^{n}$ onto a new adaptive
mesh $\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$ at time $t^{n+1}$, Section 4.2.

Step 4. Repeat Step 1 \- Step 3 until final time.

###  4.2 Adaptive Mesh Refinement/Coarsening

The purpose of mesh reconstruction is to obtain adaptive grids which delivers
higher accuracy with a lower computational cost. An efficient adaptive mesh
should have small cells in the regions with large errors and large cells in
the other regions. In particular, at time $t^{n}$, we start with the given
mesh, denoted as $\mathcal{T}^{n,m}=\\{T^{n,m}_{j}\\}$, where $T^{n,m}_{j}$ is
a triangular cell with the barycenter $(x^{n,m}_{j},y^{n,m}_{j})$, and index
$m=0,1,2...$ is the level of refinement
($\mathcal{T}^{n,0}\equiv\mathcal{T}^{0,0}$ for all $n$ and
$\mathcal{T}^{0,0}$ represents the initial mesh with no refinement). The
triangular cells in the mesh $\mathcal{T}^{n,m}$ are flagged for
refinement/de-refinement based on the weak local residual (WLR) error
estimator, see Section 4.4. On grid $\mathcal{T}^{n,m}$, we apply “regular
refinement” described in [13] on the triangles flagged for refinement to
obtain a new mesh $\mathcal{T}^{n,m+1}$ with the refinement level $m+1$.
Namely, each flagged triangle (“parent” triangle) is split into four smaller
triangles (“children” triangles) by inserting a new node at the mid-point of
each edge of the “parent” triangle. Fig. 4.1 (a) is an illustration of the
“regular refinement”, where a flagged cell $T^{n,m}_{j}$ is refined to obtain
the “children” cells $T^{n,m+1}_{j_{s}},s=1,2,3,4$ by using the mid-points of
the sides. In addition, due to the insertion of new nodes on the edges of the
non-flagged triangles adjacent to refined triangles, we must also refine these
neighboring cells by creating a new edge between the hanging node and the
opposite corner as illustrated in Fig. 4.1 (b).

(a) Triangle $T^{n,m}_{j}$ (left) is split into four “children” cells
$T^{n,m+1}_{js},s=1,2,3,4$ (right). (b) Refinement in the neighboring cells of
$T^{n,m}_{j}$.

Figure 4.1: An outline of the “regular refinement”.

In practice, there may be some cells having very large WLR error (4.9), and we
need to reach a higher level of refinement for those cells to improve the
accuracy. This can be done by repeating the refinement for the flagged
triangles in the refined mesh $\mathcal{T}^{n,m+1}$ to get the mesh with
higher level $\mathcal{T}^{n,m+2},m=0,1,2,....$. Fig. 4.2 is the illustration
of the “regular refinement” procedure with two levels of refinement.

Figure 4.2: An example of the “regular refinement” procedure with two levels
of refinement. The left figure is the initial coarse mesh $\mathcal{T}^{n,0}$
with the region flagged for the refinement (gray). The middle figure is the
“first” level mesh $\mathcal{T}^{n,1}$ with the region flagged for higher
level of the refinement (gray). The right figure is the “second” level mesh
$\mathcal{T}^{n,2}$.

Note that in the numerical simulations of the wave phenomena, the regions of
the domain flagged for refinement changes over time evolution and the
refinement in some cells may become no longer needed. Hence, in [13], the de-
refinement/coarsening procedure is performed to deactivate unnecessarily fine
cells in the grid. Namely, at time $t^{n}$, we deactivate “children” cells in
the mesh $\mathcal{T}^{n,m+1},m=0,1,...,M_{n}-1$ based on the WLR and activate
the corresponding “parent” cell from the mesh $\mathcal{T}^{n,m}$ back. More
details of refinement/de-refinement process can be seen in [13].

Finally, at time $t^{n}$, a hierarchical system of grids
$\mathcal{S}^{n}=\\{\mathcal{T}^{n,0},\mathcal{T}^{n,1},\mathcal{T}^{n,2},...,\mathcal{T}^{n,\mathcal{M}_{n}}\\}$
is obtained, where $\mathcal{T}^{n,m}$, $m=1,2,...,\mathcal{M}_{n}$ is the
grid with the level of refinement $m$ reconstructed from the grid
$\mathcal{T}^{n,m-1}$. We will use the final mesh
$\mathcal{T}^{n,\mathcal{M}_{n}}\in\mathcal{S}^{n}$ for adaptive central-
upwind scheme (4.5a)-(4.5b) to evolve the numerical solution from time $t^{n}$
to time $t^{n+1}$. Next, at $t^{n+1}$, we generate a new adaptive grid
$\mathcal{T}^{n+1,\mathcal{M}_{n+1}}\in\mathcal{S}^{n+1}$ from the mesh
$\mathcal{T}^{n,\mathcal{M}_{n}}$. After the mesh reconstruction at time
$t^{n+1}$, the obtained cell averages $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}^{n+1}$ on the mesh
$\mathcal{T}^{n,\mathcal{M}_{n}}$ need to be accurately projected on the new
mesh $\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$, using the ideas as summarized
briefly below.

Case 1. At $t^{n+1}$, a triangle
$T^{n+1,\mathcal{M}_{n+1}}_{j}\in\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$ is the
same cell as in the grid $\mathcal{T}^{n,\mathcal{M}_{n}}$, we will maintain
the cell averages for that triangle at $t^{n+1}$.

Case 2. A cell
$T^{n+1,\mathcal{M}_{n+1}}_{j}\in\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$ is
obtained by de-refining some finer cells
$T^{n,\mathcal{M}_{n}}_{j_{s}}\in\mathcal{T}^{n,\mathcal{M}_{n}},s=1,2,..,S$.
The cell average of solution, $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$U$\kern-0.50003pt}}}^{n+1}_{j}$ in the cell
$T^{n+1,\mathcal{M}_{n+1}}_{j}$, is computed as

$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$U$\kern-0.50003pt}}}^{n+1}_{j}$
$\displaystyle=\dfrac{1}{|T^{n+1,\mathcal{M}_{n+1}}_{j}|}\sum\limits_{s=1}^{S}\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$U$\kern-0.50003pt}}}^{n}_{j_{s}}|T^{n,\mathcal{M}_{n}}_{j_{s}}|,$

where $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$U$\kern-0.50003pt}}}^{n}_{j_{s}}$ is the
solution in $T^{n,\mathcal{M}_{n}}_{j_{s}}$.

Case 3. A triangle
$T^{n+1,\mathcal{M}_{n+1}}_{j}\in\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$ is
obtained from the refinement of the cell
$T^{n,\mathcal{M}_{n}}_{i}\in\mathcal{T}^{n,\mathcal{M}_{n}}$. If
$T^{n,\mathcal{M}_{n}}_{i}$ is a single-fluid cell, the cell averages at
$t^{n+1}$ in $\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$ are approximated by using
the the piecewise linear reconstruction (2.7) of the solution at $t^{n+1}$ in
the triangle $T^{n,\mathcal{M}_{n}}_{i}$. If $T^{n,\mathcal{M}_{n}}_{i}$
contains the density discontinuity, we will compute the cells averages in the
“son” cell $T^{n+1,\mathcal{M}_{n+1}}_{j}$ based on the interface tracking,
see Section 3.2, and the information of the nearby single-fluid cells. For
example, suppose that from the reconstructed interface in “dad” cell
$T^{n,\mathcal{M}_{n}}_{i}$, the “child” cell $T^{n+1,\mathcal{M}_{n+1}}_{j}$
completely lies in one fluid, called fluid 1. Hence, triangle
$T^{n+1,\mathcal{M}_{n+1}}_{j}$ is a single-fluid cell and the cell averages
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$U$\kern-0.50003pt}}}^{n+1}_{j}$ are set to be
equal to the cell averages of another “reliable” cell which is in fluid 1 and
is the closest cell to the “dad” cell $T^{n,\mathcal{M}_{n}}_{i}$.

###### Remark 4.1

The update of cell averages in case 3 does not ensure the mass conservation
which means the total mass of “son” cells in the adaptive grid is not equal to
the mass of their “dad” cell. However, the volume of fluid and the exact
location of the interface are conserved. In addition, the reconstructing
method also maintains the steady state solution for “lake at rest” situations,
see example 2 Section 5, since the values of the cell averages are obtained by
using the information from nearby “reliable” cells.

###  4.3 Second-order Adaptive Time Evolution

Note that using a global time step in the adaptive algorithm may lead to a
very small time step due to the presence of much finer cells in the mesh. To
reduce the computational cost, we consider the approach based on the adaptive
time step from [13, 11, 12, 30]. The main idea of this approach is to group
cells into different levels based on the cell sizes and applying local time
step on each level to evolve from $t^{n}$ to $t^{n+1}$. Recently, in our work
on the shallow water model [13], we have derived a simple and efficient
adaptive time evolution algorithm based on the second-order strong stability
preserving Runge-Kutta methods (SSPRK2) in [16, 11, 12, 30]. This method is
also capable to perform on the multi-fluid system (1.1a)-(1.1b). The algorithm
can be briefly described by an example as illustrated in Fig. 4.3.

Figure 4.3: The example of SSPRK2 on mesh with three cell levels, $l=0,1,2$.

First, we group all cells in the grid $\mathcal{T}^{n,\mathcal{M}_{n}}$ at
time $t^{n}$ in cell levels $l=0,1,..,L$ based on their sizes. Namely, a cell
$T^{n,\mathcal{M}_{n}}_{j}\in\mathcal{T}^{n,\mathcal{M}_{n}}$ belongs to the
level $l$, if $l$ is the smallest positive integer satisfying,

$2^{l}\geq\dfrac{\displaystyle\max_{j}\left(\displaystyle\min_{k}\left(r_{jk}\right)\right)}{\displaystyle\min_{k}(r_{j,k})},$
(4.1)

where $r_{jk},k=1,2,3$ are three altitudes of triangle
$T^{n,\mathcal{M}_{n}}_{j}$. Next, at $t^{n}$, we define the reference time
step $\Delta t$ as the local time step on the coarsest level $l=0$ of cells in
the mesh $\mathcal{T}^{n,\mathcal{M}_{n}}$ by considering the CFL-type
condition (2.20) locally on level $l=0$.

$\displaystyle\Delta t\equiv\Delta
t^{n,0}_{0}=\frac{0.9\displaystyle\max_{j}\left(\displaystyle\min_{k}\left(r_{jk}\right)\right)}{18a_{max}},$
(4.2)

where

$a_{max}:=\displaystyle\max_{j,k}(a^{\rm in}_{jk},a^{\rm out}_{jk}),$ (4.3)

and $(a^{\rm in}_{jk},a^{\rm out}_{jk})$ are the local one-sided speeds of
propagation (2.8) at $t^{n}$ for sides $k=1,2,3$ in the triangle
$T^{n,\mathcal{M}_{n}}_{j}\in\mathcal{T}^{n,\mathcal{M}_{n}}$. We set,
$t^{n+1}=t^{n}+\Delta t$.

Next, assume that $\mathcal{P}_{l}$ is the number of steps needed for the
higher levels $l=1,..L$ to evolve from $t^{n}$ to $t^{n+1},$ namely
$[t^{n},t^{n+1}]=\cup[t_{l}^{n,p},t_{l}^{n,p+1}],p=0,...,\mathcal{P}_{l}-1$
with $t^{n,0}_{l}\equiv t^{n}$, $t_{l}^{n,\mathcal{P}_{l}}\equiv
t^{n+1}\quad\forall l$. At $t^{n,p}_{l}$, the local time step for cells on
these levels $l=1,...,L$ is calculated by

$\Delta t^{n,p}_{l}=\dfrac{2^{-l}\Delta t}{\max(\mu^{n,p}_{l},1)},$ (4.4)

where parameter $\mu^{n,p}_{l}$ takes into account change in the local one-
sided speeds of the propagation,

$\mu^{n,p}_{l}=\dfrac{\displaystyle\max_{j,k}(a^{\rm in}_{jk},a^{\rm
out}_{jk})^{n,p}_{l}}{a_{max}},$

where $(a^{\rm in}_{jk},a^{\rm out}_{jk})^{n,p}_{l}$ are the local one-sided
speeds of propagation of the cell $T^{n,\mathcal{M}_{n}}_{j}$ in the level $l$
at $t^{n,p}_{l}$. Therefore, on each cell $T^{n,\mathcal{M}_{n}}_{j}$ of level
$l$, for each substep
$[t_{l}^{n,p},t_{l}^{n,p+1}]\equiv[t_{l}^{n,p},t_{l}^{n,p}+\Delta
t_{l}^{n,p}],p=0,1,2,3,...,\mathcal{P}_{l}-1$ of the evolution from $t^{n}$ to
$t^{n+1}$, we apply the following two adaptive steps of the SSPRK2 method, see
[16, 11, 12, 30],

$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{(1)}$
$\displaystyle=\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n,p}-\Delta
t^{n,p}_{l}\left(\frac{1}{|T^{n,\mathcal{M}_{n}}_{j}|}\sum_{k=1}^{3}{\bm{H}}^{n,p}_{jk}-\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{S}$\kern-0.50003pt}}}^{n,p}_{j}\right):={\bm{R}}(\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n,p},\Delta
t_{l}^{n,p}),$ (4.5a) $\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n,p+1}$
$\displaystyle=\dfrac{1}{2}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n,p}+\dfrac{1}{2}{\bm{R}}(\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{(1)},\Delta
t_{l}^{n,p}).$ (4.5b)

The flux term ${\bm{H}}^{n,p}_{jk}$ in (4.5a) -(4.5b) is the flux (2.6)
computed at $t=t^{n,p}_{l}$. The source term $\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt${\bm{S}}$\kern-0.50003pt}}}^{n,p}_{j}$ in
(4.5a) -(4.5b) is the source (2.25) - (2.26) computed at $t=t^{n,p}_{l}$ with
the time step $\Delta t_{l}^{n,p}$. Note that, $\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n,0}\equiv\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n}$ and
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n,\mathcal{P}_{l}}\equiv\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n+1}$.

###### Remark 4.2

If cells from different cell levels are neighbors, we use linear interpolation
in time to match the time levels of such cells, see also Fig. 3.6, for the
illustration of the interpolation.

###  4.4 A Posteriori Error Estimator

To create a robust indicator for the adaptive mesh refinement, in our prior
work [13], we have derived local error estimator from the idea of Weak Local
Residual (WLR) presented in [40, 18]. This error indicator has shown its
advantages in accurately capturing the regions with large error in numerical
simulation for Saint-Venant system of shallow water model. Hence, for the
adaptive central-upwind scheme for SWEDs, we will extend the error estimator
by applying the computation performed in [13] for the last equations in the
system (1.1d).

Let us recall that from the weak form of the mass conservation equation
(1.1a), in [13], the WLR errors $E^{w,n+\frac{1}{2}}_{i}$ at each node $N_{i}$
on mesh $\mathcal{T}^{n,\mathcal{M}_{n}}$ is given by the formula,

$\displaystyle E^{w,n+\frac{1}{2}}_{i}$
$\displaystyle=\frac{1}{\Delta}(\mathcal{U}^{w,n+\frac{1}{2}}_{i}+\mathcal{F}^{w,n+\frac{1}{2}}_{i}+\mathcal{G}^{w,n+\frac{1}{2}}_{i}),$
(4.6) $\displaystyle\mathcal{\bm{U}}^{w,n+\frac{1}{2}}_{i}$
$\displaystyle=\sum\limits_{c=1}^{C_{i}}\frac{1}{3}|T^{n,\mathcal{M}_{n}}_{j_{c}}|(\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$w$\kern-0.50003pt}}}^{n}_{j_{c}}-\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$w$\kern-0.50003pt}}}^{n+1}_{j_{c}}),$
$\displaystyle\mathcal{F}^{w,n+\frac{1}{2}}_{i}$
$\displaystyle=\sum\limits_{c=1}^{C_{i}}a^{(i)}_{c}\frac{\Delta
t}{2}|T^{n,\mathcal{M}_{n}}_{j_{c}}|((\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$hu$\kern-0.50003pt}}})^{n}_{j_{c}}+(\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$hu$\kern-0.50003pt}}})^{n+1}_{j_{c}}),$
$\displaystyle\mathcal{G}^{w,n+\frac{1}{2}}_{i}$
$\displaystyle=\sum\limits_{c=1}^{C_{i}}b^{(i)}_{c}\frac{\Delta
t}{2}|T^{n,\mathcal{M}_{n}}_{j_{c}}|((\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$hv$\kern-0.50003pt}}})^{n}_{j_{c}}+(\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$hv$\kern-0.50003pt}}})^{n+1}_{j_{c}}),$

where $C_{i}$ is the number of triangles $T^{n,\mathcal{M}_{n}}_{j_{c}}$
having common vertex at node $N_{i}$. Here, the quantity
$(a^{(i)}_{c},b^{(i)}_{c})$ is the gradient of the linear piece restricted to
$T^{n,\mathcal{M}_{n}}_{j_{c}}$,

$\displaystyle a^{(i)}_{c}$
$\displaystyle=\dfrac{\widetilde{y}_{2}-\widetilde{y}_{3}}{(\widetilde{y}_{3}-\widetilde{y}_{i})(\widetilde{x}_{2}-\widetilde{x}_{i})-(\widetilde{y}_{2}-\widetilde{y}_{i})(\widetilde{x}_{3}-\widetilde{x}_{i})},$
(4.7) $\displaystyle b^{(i)}_{c}$
$\displaystyle=\dfrac{\widetilde{x}_{3}-\widetilde{x}_{2}}{(\widetilde{y}_{3}-\widetilde{y}_{i})(\widetilde{x}_{2}-\widetilde{x}_{i})-(\widetilde{y}_{2}-\widetilde{y}_{i})(\widetilde{x}_{3}-\widetilde{x}_{i})},$

where
$N_{i}=(\widetilde{x}_{i},\widetilde{y}_{i}),(\widetilde{x}_{2},\widetilde{y}_{2}),$
and $(\widetilde{x}_{3},\widetilde{y}_{3})$ are the three vertices of triangle
$T^{n,\mathcal{M}_{n}}_{j_{c}}$.

Now, by applying the same calculation in [13] on the weak form of the last
equations in the system (1.1d), we then define the WLR error of variable
$h\rho$ at node $N_{i}$ in the grid as,

$\displaystyle E^{h\rho,n+\frac{1}{2}}_{i}$
$\displaystyle=\frac{1}{\Delta}(\mathcal{U}^{h\rho,n+\frac{1}{2}}_{i}+\mathcal{F}^{h\rho,n+\frac{1}{2}}_{i}+\mathcal{G}^{h\rho,n+\frac{1}{2}}_{i}),$
(4.8) $\displaystyle\mathcal{\bm{U}}^{h\rho,n+\frac{1}{2}}_{i}$
$\displaystyle=\sum\limits_{c=1}^{C_{i}}\frac{1}{3}|T^{n,\mathcal{M}_{n}}_{j_{c}}|(\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}^{n}_{j_{c}}-\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}^{n+1}_{j_{c}}),$
$\displaystyle\mathcal{F}^{h\rho,n+\frac{1}{2}}_{i}$
$\displaystyle=\sum\limits_{c=1}^{C_{i}}a^{(i)}_{c}\frac{\Delta
t}{2}|T^{n,\mathcal{M}_{n}}_{j_{c}}|((\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$hu\rho$\kern-0.50003pt}}})^{n}_{j_{c}}+(\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$hu\rho$\kern-0.50003pt}}})^{n+1}_{j_{c}}),$
$\displaystyle\mathcal{G}^{h\rho,n+\frac{1}{2}}_{i}$
$\displaystyle=\sum\limits_{c=1}^{C_{i}}b^{(i)}_{c}\frac{\Delta
t}{2}|T^{n,\mathcal{M}_{n}}_{j_{c}}|((\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$hv\rho$\kern-0.50003pt}}})^{n}_{j_{c}}+(\hbox{\vbox{\hrule
height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$hv\rho$\kern-0.50003pt}}})^{n+1}_{j_{c}}).$

Hence, the error in a cell
$T^{n,\mathcal{M}_{n}}_{j}\in\mathcal{T}^{n.\mathcal{M}_{n}}$ takes into
account both WLR errors of water surface $w$ and of variable $h\rho$ as,

$e_{j}=\max_{k}\left(\left|E^{w,n+\frac{1}{2}}_{jk}\right|,\left|E^{h\rho,n+\frac{1}{2}}_{jk}\right|\right),\quad
k=1,2,3,$ (4.9)

where $E^{w,n+\frac{1}{2}}_{jk}$ and $E^{h\rho,n+\frac{1}{2}}_{jk}$ are the
WLR errors computed in (4.6) and (4.8) at a node $k$ of triangle $T_{j}$.

The error $e_{j}$ in each cell
$T^{n.\mathcal{M}_{n}}_{j}\in\mathcal{T}^{n,\mathcal{M}_{n}}$ is compared to
the error tolerance computed by

$\omega=\sigma\max_{j}(e_{j}),$ (4.10)

where $\sigma<1$ is a given problem-dependent constant. Based on the error
comparison, the cell is then either “flagged” for refinement/de-refinement or
“no-change”.

##  5 Numerical Experiments

In this section, we will verify the computational efficiency of the designed
adaptive central-upwind scheme. We compare the results of the adaptive method
with the results of the central-upwind scheme without the adaptivity, see
Section 4, on uniform triangular meshes (example of such uniform triangular
mesh is outlined in Fig. 5.1). To this end, in all examples, we calculate the
$L^{1}$-errors and a ratio,
$\mathcal{R}_{CPU}=\frac{CPU_{uniform}}{CPU_{adaptive}}$, which is the ratio
of the CPU times of the central-upwind algorithm without the mesh
reconstruction to the CPU time of the adaptive central-upwind algorithm. To
compare $L^{1}$-errors, as well as to compare the CPU times and to compute
$\mathcal{R}_{CPU}$, we consider uniform mesh and the adaptive mesh with the
same size of the smallest cells. Namely, in Tables 5.1-5.8, $L^{1}$-errors and
$\mathcal{R}_{CPU}$ are computed by using the uniform meshes $2\times N\times
N,N=100,200,400$ and the corresponding adaptive meshes which are reconstructed
from the coarser uniform mesh $2\times N/2^{\mathcal{M}}\times
N/2^{\mathcal{M}}$ ($\mathcal{M}=1,2$ is the highest level of refinement in
the adaptive mesh). To compute the $L^{1}$-errors, the reference solution is
obtained by applying the central-upwind method without implementing adaptivity
techniques on the uniform mesh with $2\times 800\times 800$ triangles. In all
experiments, we consider a zero-order extrapolation at all boundaries. In
addition, we use the gravitational acceleration, $g=1.0$ and the reference
density, $\rho_{0}=997$ [15] for all examples. The desingularization
parameters $\tau$ and $\varepsilon$ for calculations of the velocity
components $u$ and $v$ are set
$\tau=\max_{j}\\{|T^{n,\mathcal{M}_{n}}_{j}|^{2}\\}$ and $\varepsilon=10^{-4}$
(see Section 2.1 formula (2.6) in [29]).

Figure 5.1: An outline of uniform triangular mesh.

###  5.1 Example 1:

In the first example taken from [9], we will compare the performance of the
adaptive central-upwind scheme and the central-upwind scheme without
adaptivity on uniform triangular meshes. We also verify experimentally the
advantages of the interface reconstruction in compressing the diffusion of
variable density at the interface of fluids and improving the accuracy of
computed solutions.

We consider the bottom topography $B(x,y,0):=0$ and the following initial
condition,

$(w,u,v,\rho)^{T}(x,y,0)=\begin{cases}(2,0,0,1.5\rho_{0}),\quad\mbox{if}\quad&x^{2}+y^{2}<0.5,\\\
(1,0,0,\rho_{0}),\quad&\mbox{otherwise}.\end{cases}$ (5.1)

The data is simulated in the domain $[-1,1]\times[-1,1]$. The error tolerance
(4.10) for the mesh refinement in this example is set to
$\omega=0.01\max_{j}(e_{j})$.

In Fig. 5.2 and Fig. 5.3, we show the numerical solution of the water surface
(first column) and the density (second column) at time $t=0.15$. The solution
is calculated by using the central-upwind scheme on uniform meshes on Fig. 5.2
(a, b) and Fig. 5.3 (a,b) and using the adaptive central-upwind scheme on Fig.
5.2 (c, d) and Fig. 5.3 (c,d). The adaptive meshes in Fig. 5.2 (third column)
are obtained from the uniform mesh $2\times 100\times 100$, Fig. 5.2 (a). The
adaptive mesh with one level of refinement $\mathcal{M}=1$ (as the highest
level of refinement) is on Fig. 5.2 (c), and with two levels of refinement
$\mathcal{M}=2$ (as the highest level of refinement) is on Fig. 5.2 (d). As
can be observed in Fig. 5.2, both $w$ and $\rho$ are much sharper resolved by
using the adaptive central-upwind scheme. We note also, that by increasing the
level of refinement from $\mathcal{M}=1$ to $\mathcal{M}=2$, the number of
cells in the mesh increases from $40,292$ cells with $\mathcal{M}=1$ to
$60,326$ cells with $\mathcal{M}=2$, but the accuracy is clearly improved with
higher resolution as seen in Fig. 5.2 (c) and Fig. 5.2 (d). Also from the
adaptive meshes in Fig. 5.2 (third column), one can easily see that only cells
in the region having steep gradients are refined. This means that the WLR
error estimator accurately detects regions in the domain for adaptive
refinement/coarsening.

(a) Uniform mesh $2\times 100\times 100$.

(b) Uniform mesh $2\times 200\times 200$.

(c) Adaptive mesh with $\mathcal{M}=1$.

(d) Adaptive mesh with $\mathcal{M}=2$.

Figure 5.2: Example 1: Contour plots of computational water surface
$w(x,y,0.15)$ (first column) and density $\rho(x,y,0.15)$ (second column) of
the IVP (5.1) with the corresponding meshes (third column).

(a) Uniform mesh $2\times 100\times 100$.

(b) Uniform mesh $2\times 200\times 200$.

(c) Adaptive mesh with $\mathcal{M}=1$.

(d) Adaptive mesh with $\mathcal{M}=2$.

Figure 5.3: Example 1: 3-D plots of computational water surface $w(x,y,0.15)$
(first column) and density $\rho(x,y,0.15)$ (second column) of the IVP (5.1).

Next, in Table 5.1 we calculate the $L^{1}$-errors obtained on the adaptive
grids and on the fixed uniform grids. The errors obtained in the uniform
meshes are approximate to the errors calculated in the corresponding adaptive
meshes (the adaptive meshes have the same size of the smallest cells with the
uniform meshes). However, the adaptive scheme uses fewer cells than the
central-upwind scheme which does not consider the mesh reconstruction. In
addition, in Table 5.2, we also compute the $\mathcal{R}_{CPU}$ ratio to
compare the computational cost of the two methods. The results in Table 5.1
and Table 5.2 show that the adaptive scheme produces similar accuracy as the
scheme designed on fixed uniform triangular meshes, but at a less
computational cost.

algorithm without adaptivity | adaptive algorithm
---|---
| one level $\mathcal{M}=1$ | two levels $\mathcal{M}=2$
cells |  $L^{1}$-error | rate | cells |  $L^{1}$-error | rate | cells |  $L^{1}$-error | rate
$2\times 100\times 100$ | 0.0256 |  | 13,516 | 0.0257 |  | 12,800 | 0.0265 |
$2\times 200\times 200$ | 0.0127 | 1.01 | 40,292 | 0.0128 | 1.00 | 38,284 | 0.0133 | 0.99
$2\times 400\times 400$ | 0.0047 | 1.43 | 153,152 | 0.0049 | 1.39 | 60,326 | 0.0050 | 1.41
Table 5.1: Example 1: $L^{1}$-errors of the water surface $w$ of the IVP (5.1) at $t=0.15$, and the convergence rates of the central-upwind scheme without adaptivity (uniform mesh $2\times N\times N,N=100,200,400$) and the adaptive scheme (the corresponding adaptive mesh is reconstructed from the uniform mesh $2\times N/2^{\mathcal{M}}\times N/2^{\mathcal{M}}$). uniform mesh |  adaptive mesh $\mathcal{M}=1$ |  $\mathcal{R}_{CPU}$ with $\mathcal{M}=1$ |  adaptive mesh $\mathcal{M}=2$ |  $\mathcal{R}_{CPU}$ with $\mathcal{M}=2$
---|---|---|---|---
(cells) | (cells) | total |  without grid generation | (cells) | total |  without grid generation
$2\times 100\times 100$ | 13,516 | 1.81 | 2.01 | 12,800 | 1.83 | 1.94
$2\times 200\times 200$ | 40,292 | 1.82 | 2.14 | 38,284 | 2.14 | 2.30
$2\times 400\times 400$ | 153,152 | 1.98 | 2.29 | 60,326 | 2.33 | 2.50
$\mathcal{R}_{CPU}$ average: | 1.87 | 2.15 |  | 2.10 | 2.27
Table 5.2: Example 1: $\mathcal{R}_{CPU}$ ratio for the IVP (5.1) at $t=0.15$,
where for adaptive central-upwind scheme, we consider the total CPU times and
CPU times without the grid generation.

In addition, we will use this example to show that interface reconstruction
presented in Section 3.2 plays an important part in preserving the sharpness
of the solution as well as improving the accuracy of the adaptive central-
upwind scheme. On Fig. 5.4 (a), we plot the water surface $w$ (first) and
density $\rho$ (second column) at $t=0.15$ by using the adaptive central-
upwind method, but implemented without the interface tracking technique. We
then compare the results shown in Fig. 5.4 (a) to the results calculated by
using the adaptive scheme with the interface tracking, Fig. 5.4 (b). The
adaptive meshes on Fig. 5.4 (third column) are reconstructed from the uniform
mesh $2\times 100\times 100$ with one level of refinement. As can be seen in
Fig. 5.4 (a), both $w$ and $\rho$ are very scattered around the contact wave
when we do not track and reconstruct the interface. Meanwhile, in Fig. 5.4
(b), the proposed adaptive scheme, though using an adaptive mesh with fewer
cells (10828 cells fewer), provides more accurate results.

(a) Adaptive scheme without interface tracking on adaptive mesh with 51768
cells.

(b) Adaptive scheme with interface tracking on adaptive mesh with 40940 cells

Figure 5.4: Example 1: Computational water surface $w(x,y,0.15)$ (first
column), density $\rho(x,y,0.15)$ (second column) of the IVP (5.1), and the
corresponding adaptive meshes (third column) obtained by using the proposed
adaptive central-upwind scheme (bottom) and using the adaptive scheme without
interface tracking (top).

In the next numerical test, we replace the flat bottom with the bottom
topography that consists of two Gaussian shaped humps as

$B(x,y,t)=\begin{cases}0.5e^{-100(x+0.5)^{2}+(y+0.5)^{2}},\quad&\mbox{if}\quad
x<0,\\\ 0.6e^{-100(x-0.5)^{2}+(y-0.5)^{2}},\quad&\mbox{if}\quad x\geq
0.\end{cases}$ (5.2)

The purpose of this test is to illustrate the performance of the adaptive
algorithm in situations having irregular bottom topography. In Fig. 5.5, we
show the contour plots of the water surface $w$ (first column) and density
$\rho$ (second column) obtained at $t=0.2$ by using the central-upwind scheme
with and without adaptivity. The computed solutions of the water surface
exhibit reflecting waves where the flow meets the submerged humps. Clearly,
from the plots of the adaptive meshes in Fig. 5.5 (third column), the meshes
are adapted to the behavior of the flow. Hence, the WLR error estimator is
capable to exactly detect the location of the steep local gradients in the
solution.

(a) Uniform mesh $2\times 100\times 100$.

(b) Uniform mesh $2\times 200\times 200$.

(c) Adaptive mesh with $\mathcal{M}=1$.

(d) Adaptive mesh with $\mathcal{M}=2$.

Figure 5.5: Example 1: Contour plots of computational water surface
$w(x,y,0.2)$ (first column) and density $\rho(x,y,0.2)$ (second column) of the
IVP (5.1)-(5.2) with the corresponding meshes (third column).

We then recompute the accuracy of the solution for this example 5.1-5.2, see
Table 5.3, and the CPU time ratio, see Table 5.4. The results show that the
adaptive scheme uses fewer cells and takes a smaller CPU time to achieve the
approximately small $L^{1}$-errors as computed by the scheme without
adaptivity. Therefore, the advantages of the adaptive scheme is maintained in
examples with irregular bottom level.

algorithm without adaptivity | adaptive algorithm
---|---
| one level $\mathcal{M}=1$ | two levels $\mathcal{M}=2$
cells |  $L^{1}$-error | rate | cells |  $L^{1}$-error | rate | cells |  $L^{1}$-error | rate
$2\times 100\times 100$ | 0.0271 |  | 16,252 | 0.0273 |  | 16,006 | 0.0282 |
$2\times 200\times 200$ | 0.0134 | 1.01 | 57,300 | 0.0136 | 1.00 | 42,602 | 0.0140 | 0.99
$2\times 400\times 400$ | 0.0050 | 1.43 | 213,268 | 0.0054 | 1.40 | 120,654 | 0.0055 | 1.41
Table 5.3: Example 1: $L^{1}$-errors of the water surface $w$ of the IVP (5.1)-(5.2) at $t=0.2$, and the convergence rates of the central-upwind scheme without adaptivity (uniform mesh $2\times N\times N,N=100,200,400$) and the adaptive scheme (the corresponding adaptive mesh is reconstructed from the uniform mesh $2\times N/2^{\mathcal{M}}\times N/2^{\mathcal{M}}$). uniform mesh |  adaptive mesh $\mathcal{M}=1$ |  $\mathcal{R}_{CPU}$ with $\mathcal{M}=1$ |  adaptive mesh $\mathcal{M}=2$ |  $\mathcal{R}_{CPU}$ with $\mathcal{M}=2$
---|---|---|---|---
(cells) | (cells) | total |  without grid generation | (cells) | total |  without grid generation
$2\times 100\times 100$ | 16,252 | 1.69 | 1.89 | 16,006 | 1.84 | 1.94
$2\times 200\times 200$ | 57,300 | 2.57 | 2.95 | 42,602 | 3.41 | 3.68
$2\times 400\times 400$ | 213,268 | 2.24 | 2.60 | 120,654 | 2.63 | 2.84
$\mathcal{R}_{CPU}$ average: | 2.17 | 2.48 |  | 2.63 | 2.82
Table 5.4: Example 1: $\mathcal{R}_{CPU}$ ratio for the IVP (5.1)-(5.2) at
$t=0.2$, where for adaptive central-upwind scheme, we consider the total CPU
times and CPU times without the grid generation.

###  5.2 Example 2:

The second numerical example here was proposed in [9] to verify the capability
of the adaptive scheme in preserving the steady state solution in “lake at
rest’ problems, (1.2) and (1.3). In particular, the initial data consists of
two “lake at rest” states of type (1.2) connected through the density jump
corresponding to the “lake at rest” state of type (1.3) as

$(w,u,v,\rho)^{T}(x,y,0)=\begin{cases}(3,0,0,\dfrac{4}{3}\rho_{0}),\quad\mbox{if}\quad&x^{2}+y^{2}<0.25,\\\
(2,0,0,3\rho_{0}),\quad&\mbox{otherwise}.\end{cases}$ (5.3)

In this example, we consider the bottom topography (5.2) on a computational
domain $[-1,1]\times[-1,1]$. To reconstruct the adaptive meshes, the threshold
is set, $\omega=0.1\max_{j}(e_{j})$. In Fig. 5.6, we present the plots of the
computed water surface (first column) and density (second column) at $t=0.15$
obtained by using the central-upwind scheme, but without adaptivity, Fig. 5.6
(a, b) and by using the adaptive algorithm Fig. 5.6 (c, d). The adaptive grids
plotted on Fig. 5.6 (third column) are generated from the uniform mesh
$2\times 100\times 100$ with one level of refinement $\mathcal{M}=1$, Fig. 5.6
(c), and two levels of refinement $\mathcal{M}=2$, Fig. 5.6 (d). Fig. 5.7
shows the 3D plots of the numerical solution computed by the two methods. As
expected, in Fig. 5.6 and Fig. 5.7, the adaptive scheme with interface
tracking exactly preserves the steady state. Hence, in Fig. 5.6 (third
column), the WLR error only marks cells surrounding the circle of density jump
for refinement. In addition, no pressure oscillations are observed at the
interface.

(a) Uniform mesh $2\times 100\times 100$.

(b) Uniform mesh $2\times 200\times 200$.

(c) Adaptive mesh with $\mathcal{M}=1$.

(d) Adaptive mesh with $\mathcal{M}=2$.

Figure 5.6: Example 2: Contour plots of computational water surface
$w(x,y,0.15)$ (first column) and density $\rho(x,y,0.15)$ with the
corresponding meshes (right column).

(a) Uniform mesh $2\times 100\times 100$.

(b) Uniform mesh $2\times 200\times 200$.

(c) Adaptive mesh with $\mathcal{M}=1$.

(d) Adaptive mesh with $\mathcal{M}=2$.

Figure 5.7: Example 2: 3-D plots of computational water surface $w(x,y,0.15)$
(left column) and density $\rho(x,y,0.15)$ (right column).

Next, we will illustrate the advantages of the adaptive central-upwind scheme.
We compute the $L^{1}$-error, Table 5.5, and the CPU ratio, Table 5.6, by
using the central-upwind method without adaptivity and using the adaptive
algorithm presented in our work. In order to compare the computational costs
and calculate $\mathcal{R}_{CPU}$, we consider uniform and adaptive meshes
with the same size of the smallest cells. Table 5.5 shows that in the adaptive
meshes, we achieve $L^{1}$-errors as small as the errors obtained in the
corresponding uniform meshes. However, the adaptive algorithm uses fewer cells
and reduces the CPU times up to eight times. The computational cost is
remarkably cut down since as illustrated in Fig. 5.6 and Fig. 5.7, only a few
cells in the neighborhood of the density discontinuity have large WLR errors
and are therefore marked for refinement. This example has clearly show the
efficiency of the proposed scheme for numerically solving the system of multi-
fluid flow.

algorithm without adaptivity | adaptive algorithm
---|---
| one level $\mathcal{M}=1$ | two levels $\mathcal{M}=2$
cells |  $L^{1}$-error | rate | cells |  $L^{1}$-error | rate | cells |  $L^{1}$-error | rate
$2\times 100\times 100$ | 0.0045 |  | 6,284 | 0.0045 |  | 5,928 | 0.0045 |
$2\times 200\times 200$ | 0.0021 | 1.10 | 22,546 | 0.0021 | 1.10 | 11,260 | 0.0021 | 1.10
$2\times 400\times 400$ | 7.9252e-04 | 1.41 | 85,132 | 8.5079e-04 | 1.30 | 33,904 | 8.9818e-04 | 1.23
Table 5.5: Example 2: $L^{1}$-errors of the water surface $w$ at $t=0.15$, and the convergence rates of the central-upwind scheme without adaptivity (uniform mesh $2\times N\times N,N=100,200,400$) and the adaptive scheme (the corresponding adaptive mesh is reconstructed from the uniform mesh $2\times N/2^{\mathcal{M}}\times N/2^{\mathcal{M}}$). uniform mesh |  adaptive mesh $\mathcal{M}=1$ |  $\mathcal{R}_{CPU}$ with $\mathcal{M}=1$ |  adaptive mesh $\mathcal{M}=2$ |  $\mathcal{R}_{CPU}$ with $\mathcal{M}=2$
---|---|---|---|---
(cells) | (cells) | total |  without grid generation | (cells) | total |  without grid generation
$2\times 100\times 100$ | 6,284 | 3.01 | 3.54 | 5,928 | 3.48 | 3.75
$2\times 200\times 200$ | 22,546 | 3.78 | 4.76 | 11,260 | 7.76 | 8.66
$2\times 400\times 400$ | 85,132 | 4.46 | 5.57 | 33,904 | 11.13 | 12.86
$\mathcal{R}_{CPU}$ average: | 3.75 | 4.62 |  | 7.46 | 8.42
Table 5.6: Example 2: $\mathcal{R}_{CPU}$ ratio at $t=0.15$, where for
adaptive central-upwind scheme, we consider the total CPU times and CPU times
without the grid generation.

###  5.3 Example 3:

The last example is designed to illustrate the capability of the proposed
adaptive algorithm to handle irregular density interfaces. Hence, in a domain
$[-1,1]\times[-1,1]$, the density jump at $t=0$ is given by a curve which
consists of a horizontal segment, a vertical segment, and a quarter of a
circle connected at their endpoints as illustrated in Fig. 5.8.

Figure 5.8: Example 3: The interface at initial time $t=0$.

The initial condition is described by

$(w,u,v,\rho)^{T}(x,y,0)=\begin{cases}(2,0,0,\rho_{0}),\quad\mbox{if}\quad&(x,y)\in\Omega,\\\
(1,0,0,1.5\rho_{0}),\quad&\mbox{otherwise}.\end{cases}$ (5.4)

where

$\Omega:=\\{x<-0.5,y<0\\}\cup\\{(x+0.5)^{2}+(y+0.5)^{2}<0.25\\}\cup\\{x<0,y<-0.5\\}.$

. We consider a bottom topography with a surbmerged hump as

$B(x,y,t)=0.5e^{-100(x^{2}+y^{2})}.$

In this example, we will also perform the same numerical tests which are done
in previous examples. Namely, we first calculate the water surface and density
at $t=0.15$ using central-upwind scheme without adaptivity and present the
results in Fig. 5.9 (a, b) and Fig. 5.10 (a, b). In 5.9 (c, d) and 5.10 (c,
d), we plot the results for $w$ (first column) and $\rho$ (second column)
obtained by the adaptive scheme. The adaptive grids in 5.9 (third column) are
generated from the uniform grid $2\times 100\times 100$ for one level of
refinement $\mathcal{M}=1$ and $\mathcal{M}=2$ using the threshold
$\omega=0.01\max_{j}(e_{j})$. As expected, the solutions obtained in the
adaptive meshes with high levels of refinement are much sharper than the
solutions computed in fixed uniform meshes. The density jump moves Northeast
and does not diffuse. There is no non-physical spurious waves generated at the
interface. Also, as can be seen in the adaptive meshes, the WLR error
indicator captures subtle features of the solution.

(a) Uniform mesh $2\times 100\times 100$.

(b) Uniform mesh $2\times 200\times 200$.

(c) Adaptive mesh with $\mathcal{M}=1$.

(d) Adaptive mesh with $\mathcal{M}=2$.

Figure 5.9: Example 3: Contour plots of computational water surface
$w(x,y,0.15)$ (first column) and density $\rho(x,y,0.15)$ (second column) with
the corresponding meshes (third column).

(a) Uniform mesh $2\times 100\times 100$.

(b) Uniform mesh $2\times 200\times 200$.

(c) Adaptive mesh with $\mathcal{M}=1$.

(d) Adaptive mesh with $\mathcal{M}=2$.

Figure 5.10: Example 3: 3-D plots of computational water surface $w(x,y,0.15)$
(left column) and density $\rho(x,y,0.15)$ (right column).

Finally, we compute the $L^{1}$-errors in Table 5.7 and the CPU ratio in Table
5.8 by using the adaptive scheme and using the central-upwind method without
the mesh generation. By comparing the results presented in Tables 5.7 and 5.8,
one can easily see that at a reduced computational cost, the proposed adaptive
central-upwind method is still able to obtain the accurate solutions for this
example. In our experiments, we considered only $\mathcal{M}=1$ and
$\mathcal{M}=2$, but to further enhance the accuracy of the numerical solution
at the lower computational cost, one can consider higher levels of refinement.

algorithm without adaptivity | adaptive algorithm
---|---
| one level $\mathcal{M}=1$ | two levels $\mathcal{M}=2$
cells |  $L^{1}$-error | rate | cells |  $L^{1}$-error | rate | cells |  $L^{1}$-error | rate
$2\times 100\times 100$ | 0.0155 |  | 12,164 | 0.0145 |  | 8,994 | 0.0166 |
$2\times 200\times 200$ | 0.0074 | 1.07 | 40,814 | 0.0074 | 0.97 | 30,263 | 0.0071 | 1.23
$2\times 400\times 400$ | 0.0027 | 1.45 | 148,473 | 0.0028 | 1.40 | 87,214 | 0.0028 | 1.34
Table 5.7: Example 3: $L^{1}$-errors of the water surface $w$ at $t=0.15$, and the convergence rates of the central-upwind scheme without adaptivity (uniform mesh $2\times N\times N,N=100,200,400$) and the adaptive scheme (the corresponding adaptive mesh is reconstructed from the uniform mesh $2\times N/2^{\mathcal{M}}\times N/2^{\mathcal{M}}$). uniform mesh |  adaptive mesh $\mathcal{M}=1$ |  $\mathcal{R}_{CPU}$ with $\mathcal{M}=1$ |  adaptive mesh $\mathcal{M}=2$ |  $\mathcal{R}_{CPU}$ with $\mathcal{M}=2$
---|---|---|---|---
(cells) | (cells) | total |  without grid generation | (cells) | total |  without grid generation
$2\times 100\times 100$ | 12,164 | 1.75 | 1.98 | 8,994 | 2.69 | 2.87
$2\times 200\times 200$ | 40,814 | 2.71 | 3.14 | 30,263 | 3.57 | 3.86
$2\times 400\times 400$ | 148,473 | 2.86 | 3.38 | 87,214 | 4.60 | 5.03
$\mathcal{R}_{CPU}$ average: | 2.44 | 2.83 |  | 3.62 | 3.92
Table 5.8: Example 3: $\mathcal{R}_{CPU}$ ratio at $t=0.15$, where for
adaptive central-upwind scheme, we consider the total CPU times and CPU times
without the grid generation.

##  6 Conclusion

We have developed a new adaptive well-balanced and positivity preserving
central-upwind scheme on unstructured traingular meshes for shallow water
equations with variable density. The scheme is designed as an extension the
scheme in [13] by utilizing the interface tracking method in [9] and the
interface reconstruction in [15]. The proposed scheme is capable to preserve
the steady state solutions (1.2) and (1.3) and prevent the oscillation at the
density jumps. In addition, to achieve an efficient strategy for the adaptive
mesh reconstruction, we also obtain a robust local error indicator. We
performed several challenging numerical tests for multi-fluid models and we
demonstrated that the new adaptive central-upwind scheme maintains well-
balanced and positivity-preserving properties and obtains high-accuracy at a
reduced computational cost.

## Acknowledgements

I would like to sincerely thank my advisor, Prof. Yekaterina Epshteyn, for her
suggestion of the problem to investigate, her encouragement and help. The
proposed approach in this paper is an extension of the adaptive method
developed in [13] which is my joint work with her. Without her support, this
work would not have started, progressed, or ended.

## References

##

  * [1] D. S. Balsara, M. Dumbser, and R. Abgrall, Multidimensional hllc riemann solver for unstructured meshes–with application to euler and mhd flows, Journal of Computational Physics, 261 (2014), pp. 172–208.
  * [2] A. Bollermann, G. Chen, A. Kurganov, and S. Noelle, A well-balanced reconstruction of wet/dry fronts for the shallow water equations, J. Sci. Comput., 56 (2013), pp. 267–290.
  * [3] S. Bryson, Y. Epshteyn, A. Kurganov, and G. Petrova, Central Upwind Scheme on Triangular Grids for the Saint Venant System of Shallow Water Equations, AIP Conference Proceedings, 1389 (2011), pp. 686–689.
  * [4] S. Bryson, Y. Epshteyn, A. Kurganov, and G. Petrova, Well-balanced positivity preserving central-upwind scheme on triangular grids for the Saint-Venant system, ESAIM Math. Model. Numer. Anal., 45 (2011), pp. 423–446.
  * [5] S. Bryson and D. Levy, Balanced central schemes for the shallow water equations on unstructured grids, SIAM J. Sci. Comput., 27 (2005), pp. 532–552 (electronic).
  * [6] A. Chertock, Y. Epshteyn, H. Hu, and A. Kurganov, High-order positivity-preserving hybrid finite-volume-finite-difference methods for chemotaxis systems, Adv. Comput. Math., 44 (2018), pp. 327–350.
  * [7] A. Chertock, K. Fellner, A. Kurganov, A. Lorz, and P. A. Markowich, Sinking, merging and stationary plumes in a coupled chemotaxis-fluid model: a high-resolution numerical approach, Journal of Fluid Mechanics, 694 (2012), pp. 155–190.
  * [8] A. Chertock, A. Kurganov, and Y. Liu, Central-upwind schemes for the system of shallow water equations with horizontal temperature gradients, Numer. Math., 127 (2014), pp. 595–639.
  * [9] , Central-upwind schemes for the system of shallow water equations with horizontal temperature gradients, Numerische Mathematik, 127 (2014), pp. 595–639.
  * [10] P. J. Dellar, Common hamiltonian structure of the shallow water equations with horizontal temperature gradients and magnetic fields, Physics of Fluids, 15 (2003), pp. 292–297.
  * [11] M. O. Domingues, S. M. Gomes, O. Roussel, and K. Schneider, An adaptive multiresolution scheme with local time stepping for evolutionary pdes, Journal of Computational Physics, 227 (2008), pp. 3758–3780.
  * [12] R. Donat, M. Martí, A. Martínez-Gavara, and P. Mulet, Well-balanced adaptive mesh refinement for shallow water flows, Journal of Computational Physics, 257 (2014), pp. 937–953.
  * [13] Y. Epshteyn and T. Nguyen, Adaptive central-upwind scheme on triangular grids for the saint-venant system, arXiv preprint arXiv:2011.06143, (2020).
  * [14] M. A. Ghazizadeh and A. Mohammadian, An adaptive central-upwind scheme on quadtree grids for variable density shallow water equations, arXiv preprint arXiv:2008.02111, (2020).
  * [15] M. A. Ghazizadeh, A. Mohammadian, and A. Kurganov, An adaptive well-balanced positivity preserving central-upwind scheme on quadtree grids for shallow water equations, Computers $\&$ Fluids, 208 (2020), p. 104633.
  * [16] S. Gottlieb, C.-W. Shu, and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev., 43 (2001), pp. 89–112.
  * [17] C. W. Hirt and B. D. Nichols, Volume of fluid (vof) method for the dynamics of free boundaries, Journal of computational physics, 39 (1981), pp. 201–225.
  * [18] S. Karni and A. Kurganov, Local error analysis for approximate solutions of hyperbolic conservation laws, Adv. Comput. Math., 22 (2005), pp. 79–99.
  * [19] A. Kurganov, Finite-volume schemes for shallow-water equations, Acta Numerica, 27 (2018), pp. 289–351.
  * [20] A. Kurganov and D. Levy, Central-upwind schemes for the Saint-Venant system, M2AN Math. Model. Numer. Anal., 36 (2002), pp. 397–425.
  * [21] A. Kurganov and C.-T. Lin, On the reduction of numerical dissipation in central-upwind schemes, Commun. Comput. Phys, 2 (2007), pp. 141–163.
  * [22] A. Kurganov and Y. Liu, New adaptive artificial viscosity method for hyperbolic systems of conservation laws, J. Comput. Phys., 231 (2012), pp. 8114–8132.
  * [23] A. Kurganov and J. Miller, Central-upwind scheme for Savage-Hutter type model of submarine landslides and generated tsunami waves, Comput. Methods Appl. Math., 14 (2014), pp. 177–201.
  * [24] A. Kurganov, S. Noelle, and G. Petrova, Semi-discrete central-upwind schemes for hyperbolic conservation laws and Hamilton-Jacobi equations, SIAM J. Sci. Comput., 23 (2001), pp. 707–740.
  * [25] A. Kurganov, S. Noelle, and G. Petrova, Semidiscrete central-upwind schemes for hyperbolic conservation laws and hamilton–jacobi equations, SIAM Journal on Scientific Computing, 23 (2001), pp. 707–740.
  * [26] A. Kurganov and G. Petrova, Central-upwind schemes on triangular grids for hyperbolic systems of conservation laws, Numer. Methods Partial Differential Equations, 21 (2005), pp. 536–552.
  * [27] A. Kurganov and G. Petrova, A second-order well-balanced positivity preserving scheme for the Saint-Venant system, Commun. Math. Sci., 5 (2007), pp. 133–160.
  * [28] A. Kurganov and E. Tadmor, New high-resolution semi-discrete central schemes for Hamilton-Jacobi equations, J. Comput. Phys., 160 (2000), pp. 720–742.
  * [29] X. Liu, J. Albright, Y. Epshteyn, and A. Kurganov, Well-balanced positivity preserving central-upwind scheme with a novel wet/dry reconstruction on triangular grids for the Saint-Venant system, Journal of Computational Physics, 374 (2018), pp. 213 – 236.
  * [30] M. Moreira Lopes, M. O. Domingues, K. Schneider, and O. Mendes, Local time-stepping for adaptive multiresolution using natural extension of runge–kutta methods, Journal of Computational Physics, 382 (2019), pp. 291 – 318.
  * [31] S. Osher and R. P. Fedkiw, Level set methods: an overview and some recent results, Journal of Computational physics, 169 (2001), pp. 463–502.
  * [32] T. Qin, C.-W. Shu, and Y. Yang, Bound-preserving discontinuous galerkin methods for relativistic hydrodynamics, Journal of Computational Physics, 315 (2016), pp. 323–347.
  * [33] W. Rider and D. Kothe, Stretching and tearing interface tracking methods, in 12th computational fluid dynamics conference, 1995, p. 1717.
  * [34] W. J. Rider and D. B. Kothe, Reconstructing volume tracking, Journal of computational physics, 141 (1998), pp. 112–152.
  * [35] P. Ripa, Conservation laws for primitive equations models with inhomogeneous layers, Geophysical & Astrophysical Fluid Dynamics, 70 (1993), pp. 85–111.
  * [36] , On improving a one-layer ocean model with thermodynamics, Journal of Fluid Mechanics, 303 (1995), pp. 169–201.
  * [37] M. L. Sæ tra, A. R. Brodtkorb, and K.-A. Lie, Efficient GPU-implementation of adaptive mesh refinement for the shallow-water equations, J. Sci. Comput., 63 (2015), pp. 23–48.
  * [38] J. A. Sethian and P. Smereka, Level set methods for fluid interfaces, Annual review of fluid mechanics, 35 (2003), pp. 341–372.
  * [39] H. Shirkhani, A. Mohammadian, O. Seidou, and A. Kurganov, A well-balanced positivity-preserving central-upwind scheme for shallow water equations on unstructured quadrilateral grids, Comput. & Fluids, 126 (2016), pp. 25–40.
  * [40] E. Tadmor, Local error estimates for discontinuous solutions of nonlinear hyperbolic equations, SIAM J. Numer. Anal., 28 (1991), pp. 891–906.
  * [41] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, and Y.-J. Jan, A front-tracking method for the computations of multiphase flow, Journal of computational physics, 169 (2001), pp. 708–759.
  * [42] S. O. Unverdi and G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, Journal of computational physics, 100 (1992), pp. 25–37.
  * [43] Y. Xing and X. Zhang, Positivity-preserving well-balanced discontinuous galerkin methods for the shallow water equations on unstructured triangular meshes, Journal of Scientific Computing, 57 (2013), pp. 19–41.
  * [44] X. Yang and A. J. James, Analytic relations for reconstructing piecewise linear interfaces in triangular and tetrahedral grids, Journal of computational physics, 214 (2006), pp. 41–54.
  * [45] X. Yang, A. J. James, J. Lowengrub, X. Zheng, and V. Cristini, An adaptive coupled level-set/volume-of-fluid interface capturing method for unstructured triangular grids, Journal of Computational Physics, 217 (2006), pp. 364–394.