Entropy-satisfying scheme for a hierarchy of
dispersive reduced models of free surface flow

Martin Parisot*¹

¹INRIA, CNRS, Sorbonne Université, Université Paris-Diderot SPC,
Laboratoire Jacques-Louis Lions, project-team ANGE, F-75012 Paris,
France

March 10, 2021

Abstract

This work is devoted to the numerical resolution in the multidimensional
framework of a hierarchy of reduced models of the water wave equations, such
as the Serre-Green-Naghdi model. A particular attention is paid to the dissi-
pation of mechanical energy at the discrete level, that act as a stability argu-
ment of the scheme, even with source terms such space and time variation of
the bathymetry. In addition, the analysis leads to a natural way to deal with dry
areas without leakage of energy. To illustrate the accuracy and the robustness
of the strategy, several numerical experiments are carried out. In particular, the
strategy is capable of treating dry areas without special treatment.

**1 Introduction**

The propagation of surface waves is an essential issue for many applications such as port planning, tsunami propagation or marine energies. The dynamics of an incompressible, homogeneous and inviscid fluid is governed by the monovalued free surface Euler model (E), also known as water wave equations. Unfortunately, (E) is too complex to be simulated at the scale of applications and some processes such as the wet/dry front or the breaking wave can be numerically unstable. For geophysical applications, reduced models are widely used because they are easier and faster to solve numerically. The non-linear Shallow Water model (SW) is the simplest and most commonly used model nowadays. However, it is well known that (SW) is not a satisfactory model for wave propagation since dispersive effects are neglected. More precisely, the main assumptions for deriving (SW) from (E) are the so-called hydrostatic pressure in the fluid $H_{yp}^p$) and the homogeneity of the horizontal velocity in

<footnote>*martin.parisot@inria.fr</footnote>
the water column $H_{yp}^u)$. These assumptions are described in more detail in §2.2.

This work is devoted to the numerical resolution of the models that bypass
$H_{yp}^p$), i.e. takes hydrodynamic pressure into account. In the past, some reduced
models have been derived from (E) to circumvent this assumption. The Serre-Green-
Naghdi model (GN), see [22, 32, 35] is probably the most famous. Several simplified
models of it are proposed in the literature and in the following we also consider the
non-hydrostatic model (NH) [6] since it is easily classifiable in a hierarchy explained
in §2.2. All these models have the particularity of being dispersive, even without
source term, and we call them dispersive models in the following. Several works in
the literature deal with the numerical resolution of these models. In [4, 9, 10], a nu-
merical scheme is proposed based on a compact form of equations, i.e. where the
only unknowns are the water depth and the horizontal velocity, using a splitting be-
tween the shallow water equations (SW) and the dispersive part of the equations.
Despite the splitting strategy, third order derivatives still required to be discretized.
In [1], the authors propose a developed form of (NH), where the vertical-averaged
vertical velocity is introduced. This new unknown added an equation to the hy-
perbolic step and limited the dispersive step to a second-order elliptical pressure-
based equation. In the current work, the elliptical equation is solved by a velocity-
based scheme because it seems more relevant to deal with dry areas. Other strate-
gies based on a new set of variables was proposed in the literature. For example in
[27] Lagrangian variables are used. In [16], a relaxation of the constrains are used
to obtain a hyperbolic system. However, the limit of the relaxation parameters to
obtain the original Green-Naghdi model is not clear in the numerical framework. In
addition, the entropic stability is not establish for these numerical scheme.

All the dispersive models satisfy the conservation of the mechanical energy. This
property is fundamental from the point of view of mathematics because it is an ar-
gument of stability for long time solution and for applications, particularly in renew-
able energies.

The main novelty of the paper lies in the proof of entropy stability at the discrete
level. Other models are studied in the literature to simulate two-phase flows, for
example [18] and entropy-satisfying numerical scheme was proposed in such cases
[11]. These models are generally hyperbolic and the numerical strategy for hyper-
bolic models are well-understand nowadays [5, 7]. In addition, level set strategy can
also be used in this context [36]. However, these models also account the mass con-
servation of air and are not well adapted to the case of geophysical flows. In the
case of vertical-integrated models, some works propose stability property or even
convergence results for the dispersive Korteweg-de Vries equation, i.e. for a scalaire
dispersive equation in 1D [12, 15]. To our knowledge, this work is the first numerical
scheme for dispersive systems in multi-dimensional framework where the dissipa-
tion of numerical mechanical energy is proven.

This document is organized as follows. First, the free surface Euler model (E) is briefly presented in §2 as well as the hierarchy of reduced models. Then, §3 is devoted to the description of the numerical resolution within the multidimensional
framework of each model in the hierarchy. The numerical scheme is built step by step with the complexity of the models. In §3.4, a simple treatment of the boundary conditions is presented. Eventually, several numerical illustrations in a one-dimensional frame are presented in §4 which validate the accuracy and robustness of the method.

# 2 The hierarchy of free surface flow models

## 2.1 Incompressible free surface Euler equations (E)

We consider an incompressible, homogeneous, inviscid free surface flow over a non-flat bathymetry subjected to the force of gravity. The frame $(t, x, z)$ is such that $t \ge 0$ is the time, $x \in \Omega \subset \mathbb{R}^d$ is the coordinate of the horizontal plane with $d \in \{1, 2\}$, and $z \in \mathbb{R}$ is the vertical coordinate oriented in such a way that $z$ increases upwards. The flow is supposed to satisfy the free surface incompressible Euler equations with a hypothesis of free surface monovalued. More precisely, the fluid is assumed to be contained between a given lower surface $z = B(t, x)$ design as the bottom and an unknown free surface level $z = \eta(t, x) \ge B(t, x)$, and satisfies for any $B < z \le \eta$

$$ (E) \qquad \begin{cases} \nabla \cdot u + \partial_z w = 0 \\ \partial_t u + (u \cdot \nabla) u + w \partial_z u = -\nabla p \\ \partial_t w + (u \cdot \nabla) w + w \partial_z w = -\partial_z p - g \end{cases} $$

with $\nabla$ and $\nabla \cdot$ are respectively the gradient and the divergence according to the $x$-variable, the horizontal velocity $u(t, x, z) \in \mathbb{R}^d$, the vertical velocity $w(t, x, z) \in \mathbb{R}$, the pressure in the fluid $p(t, x, z) \in \mathbb{R}$ and the gravity acceleration $g \in \mathbb{R}$. At bottom, the condition of non-penetration is taken into account, i.e.

$$ \partial_t B + u_{|_{z=B}} \cdot \nabla B = w_{|_{z=B}} $$

and at free surface, the pressure is given by $p(t, x, \eta(t, x)) = P(t, x)$ and the kinematic condition is taken into account, i. e.

$$ \partial_t \eta + u_{|_{z=\eta}} \cdot \nabla \eta = w_{|_{z=\eta}} $$

The system of equation must be completed with an initial condition $\eta(0, x) = \eta^0(x)$, $u(0, x, z) = u^0(x, z)$ and $w(0, x, z) = w^0(x, z)$ which has to satisfied the compatibility conditions with the divergence free condition and the boundary condition at bottom, i.e.

$$ (1) \qquad \nabla \cdot u^0 + \partial_z w^0 = 0 \quad \text{and} \quad \partial_t B_{|_{t=0}} + u^0_{|_{z=B}} \cdot \nabla B = w^0_{|_{z=B}} $$

For readability, we introduce the homogeneous potential

$$ \phi(t, x) = gB(t, x) + P(t, x) $$

and the hydrodynamic pressure, also called non-hydrostatic pressure, defined by

$$ q(t, x, z) = p(t, x, z) - P(t, x) - g(\eta(t, x) - z). $$

Let us recall the main physical properties of the Euler model (E).
<figure><img src="image_1.png" /><figcaption>Figure 1: Illustration of the notations. (left) Interpretation of the unknowns in the vertical plan. (right) Finit volume discretization in the horizontal plan.</figcaption></figure>

**Proposition 1.** *Assuming $\eta^0(x) - B(0, x) \ge 0$, then the solution of (E) satisfies*

1.i) the water depth positivity, i.e. $\eta(t, x) - B(t, x) \ge 0$.

1.ii) the mechanical energy conservation for a sufficiently smooth solution, i.e.

$$ \partial_t (\mathcal{E}(\eta - B) + \mathcal{H}^u + \mathcal{H}^w) + \nabla \cdot (\mathcal{G}^\eta + \mathcal{G}^u + \mathcal{G}^w + \mathcal{G}^q) = (\eta - B) \partial_t \phi + q_{|z=B} \partial_t B $$

with the potential energy

$$ (2) \qquad \mathcal{E}(h) = h \left( \phi + g \frac{h}{2} \right), $$

the kinetic energy $\mathcal{H}^\Psi = \mathcal{H}(\eta, \psi)$ with

$$ \mathcal{H}(\eta, \psi) = \int_B \frac{\eta}{2} |\psi|^2 \, dz, $$

and the flux are defined by $\mathcal{G}^\eta = \mathcal{G}(\eta, u, P + g\eta)$, $\mathcal{G}^u = \mathcal{G}(\eta, u, \frac{1}{2}\|u\|^2)$, $\mathcal{G}^w = \mathcal{G}(\eta, u, \frac{1}{2}|w|^2)$ and $\mathcal{G}^q = \mathcal{G}(\eta, u, q)$ with

$$ \mathcal{G}(\eta, u, \psi) = \int_B \eta u \, dz. $$

## 2.2 Brief description of the vertical-integrated models

In the current section, we present the hierarchy of vertical-integrated models of (E) that are considered in the sequel.

### 2.2.1 Shallow water model (SW)

First of all, we present (SW) [13] where the hydrodynamic pressure is neglected. The unknowns of (SW) are the water depth $h^{\text{SW}}(t, x) \in \mathbb{R}_+$ and the vertical-averaged
horizontal velocity $\bar{u}^{\text{SW}}(t, x) \in \mathbb{R}^d$, see Figure 1. More precisely, the unknowns of (E) are approximated by

$$\eta \approx B + h^{\text{SW}}, \quad u \approx \bar{u}^{\text{SW}}, \quad w \approx 0 \quad \text{and} \quad q \approx 0.$$

(SW) can be obtained from (E) mainly assuming that:

$H_{yp}^p$ the pressure is hydrostatic, i.e. $p(t, x, z) \approx P(t, x) + g(\eta(t, x) - z)$.

$H_{yp}^u$ the horizontal velocity is close to its vertical-average, i.e. $u(t, x, z) \approx \bar{u}(t, x)$.

A fine derivation from the Navier-Stokes equations is carried out in [19]. It is shown that (SW) is an approximation of a sufficiently smooth solution of Navier-Stokes equations in $O(\epsilon)$ where $\epsilon$ is the ratio of the vertical characteristic length to the horizontal characteristic length. A higher order approximation can be derived, but it introduces dissipative terms.

### 2.2.2 Non-hydrostatic model (NH)

In the simplest dispersive model subsequently considered (NH) [6], the vertical velocity is approached by its vertical-averaged and leads to a linear approximation of hydrodynamic pressure. The unknowns of (NH) are the water depth $h^{\text{NH}}(t, x) \in \mathbb{R}_+$, the vertical-averaged horizontal velocity $\bar{u}^{\text{NH}}(t, x) \in \mathbb{R}^d$, the vertical-averaged vertical velocity $\bar{w}^{\text{NH}}(t, x) \in \mathbb{R}$ and the hydrodynamic pressure at bottom $q_B^{\text{NH}}(t, x) \in \mathbb{R}$, see Figure 1. More precisely, the unknowns $(u, w, q)$ of (E) are approximated in the vertical direction in $\mathbb{P}_0 \times \mathbb{P}_0 \times \mathbb{P}_1$ such that

$$\eta \approx B + h^{\text{NH}}, \quad u \approx \bar{u}^{\text{NH}}, \quad w \approx \bar{w}^{\text{NH}} \quad \text{and} \quad q \approx \frac{B+h-z}{h} q_B^{\text{NH}}.$$

A fine derivation of (NH) from (E) is realized in [6].

### 2.2.3 Serre-Green-Naghdi model (GN)

To satisfy the divergence free condition, the vertical velocity approximation must be of one order higher than the horizontal velocity. In the simple case, the horizontal velocity is approached in the vertical direction by means of a constant function, then the vertical velocity by means of a linear function and the hydrodynamic pressure by means of a quadratic function. It yields to (GN) [22, 32, 35]. The unknowns of (GN) are the water depth $h^{\text{GN}}(t, x) \in \mathbb{R}_+$, the vertical-averaged horizontal velocity $\bar{u}^{\text{GN}}(t, x) \in \mathbb{R}^d$, the vertical-averaged vertical velocity $\bar{w}^{\text{GN}}(t, x) \in \mathbb{R}$, the oriented standard deviation of the vertical velocity $\sigma^{\text{GN}}(t, x) \in \mathbb{R}$, the hydrodynamic pressure at bottom $q_B^{\text{GN}}(t, x) \in \mathbb{R}$ and the vertical-averaged hydrodynamic pressure $\bar{q}^{\text{GN}}(t, x) \in \mathbb{R}$, see Figure 1. More precisely, the unknowns $(u, w, q)$ of (E) are approximated in the vertical direction in $\mathbb{P}_0 \times \mathbb{P}_1 \times \mathbb{P}_2$ such that

$$\eta \approx B + h^{\text{GN}}, \quad u \approx \bar{u}^{\text{GN}}, \quad w \approx \bar{w}^{\text{GN}} + 2\sqrt{3}\frac{z - B + \frac{h}{2}}{h}\sigma^{\text{GN}} \\ \text{and} \quad q \approx \frac{B + h - z}{h}q_B^{\text{GN}} + 3\frac{(z - B)(B + h - z)}{h^2}(2\bar{q}^{\text{GN}} - q_B^{\text{GN}}).$$
Other reconstructions of (E) unknowns from that of the Serre-Green-Naghdi model can be done in particular to satisfy a curl free condition of the divergence free condition. A finer derivation of (GN) from (E) and a more detailed description is performed in [25]. It is shown that (GN) is an approximation of (E) in the order of $O(\epsilon^2)$ and $O(a)$ with a is the ratio between the characteristic wave amplitude and the characteristic water depth.

Other dispersive models, see [41], can be also considered and treated using a similar numerical strategy presented in the following. We limit the presentation to (NH) and (GN) for readability. In addition, the exponents SW, NH and GN are removed when the model used is not confusing. More precisely, each of the following subsections is devoted to a single model and the unknowns refer to it.

For each model, an entropy-satisfying scheme based on Finite volume discretization is proposed. Let us introduce some of the notations used in the following.

A tesselation $\mathbb{T}$ of the horizontal domain $\Omega \subset \mathbb{R}^d$ is considered, composed of $N = \text{Card}(\mathbb{T})$ star-shaped control volumes, see Figure 1. One denotes by $k \in \mathbb{T}$ a cell, $\mathbb{F}_k$ the set of its faces and $|k|$ its surface area. Furthermore, for a face $f$, its length is noted $|f|$ and the neighbor of $k$ through $f$ is denoted by $\underline{k}_f$ such as $k \cup \underline{k}_f = f$. The unit normal to $f$ outward to $k$ is noted by $\mathcal{N}_k^{\underline{k}_f}$. To describe the schemes, the centered discrete operator $\nabla_k^\delta : (\mathbb{R}^N)^N \to \mathbb{R}^d$ and $\nabla_k^\delta : (\mathbb{R}^d)^N \to \mathbb{R}$ are used respectively defined by

$$ (3) \qquad \begin{aligned} \nabla_k^\delta \psi &= \frac{1}{|k|} \sum_{f \in \mathbb{F}_k} \frac{\psi_{\underline{k}_f} + \psi_k}{2} \mathcal{N}_k^{\underline{k}_f} |f| \\ \text{and} \quad \nabla_k^\delta \cdot \psi &= \frac{1}{|k|} \sum_{f \in \mathbb{F}_k} \frac{\psi_k + \psi_{\underline{k}_f}}{2} \mathcal{N}_k^{\underline{k}_f} |f|. \end{aligned} $$

These discrete operators are clearly consistent with $\nabla$ and $\nabla \cdot$ respectively and are second order. The time is discretized using an adaptive time step, i.e. $t^{n+1} = t^n + \delta_t^n$, with $\delta_t^n$ satisfying a CFL condition described below. The discrete unknowns $\psi_k^n$ is the approximation of the mean value of the unknowns $\psi$ in the cell $k$ at time $t^n$. For readability, we introduce the discrete time derivative of a variable $\psi$ by

$$ (4) \qquad \partial_t^{n+1} \psi = \frac{\psi^{n+1} - \psi^n}{\delta_t^n}. $$
# 3 Numerical schemes for the reduced models

## 3.1 Shallow water model (SW)

### 3.1.1 Governing equations of (SW) and main properties

Let us begin the description of the strategy with the well-known shallow water model (SW). It reads

$$ (SW) \qquad \partial_t U^{\text{SW}} + \nabla \cdot F^{\text{SW}}(U^{\text{SW}}) = S^{\text{SW}}(U^{\text{SW}}) \quad \text{with} \quad U^{\text{SW}} = \begin{pmatrix} h \\ h\bar{u} \end{pmatrix} \\ F^{\text{SW}}(U^{\text{SW}}) = \begin{pmatrix} h\bar{u} \\ h\bar{u} \otimes \bar{u} + \frac{g}{2}h^2 \mathbf{I}_d \end{pmatrix} \quad \text{and} \quad S^{\text{SW}}(U^{\text{SW}}) = \begin{pmatrix} 0 \\ -h\nabla\phi \end{pmatrix}. $$

The initial condition of (SW) is given by $h(0, x) = h^0(x)$ and $\bar{u}(0, x) = \bar{u}^0(x)$. Let us recall the main physical properties of (SW).

**Proposition 2.** Assuming $h^0(x) \ge 0$, then the solution of (SW) satisfies

2.i) the water depth positivity, i.e. $h(t, x) \ge 0$.
It is the counterpart of Proposition 1.i).

2.ii) the mechanical energy conservation for a sufficiently smooth solution, i.e.

$$ \partial_t (\mathcal{E} + \overline{\mathcal{K}^u}) + \nabla \cdot (\overline{\mathcal{G}}^h + \overline{\mathcal{G}^u}) = h \partial_t \phi $$

with the potential energy $\mathcal{E}(h)$ defined by (2), the kinetic energy owing from the horizontal velocity $\overline{\mathcal{K}^u} = \overline{\mathcal{K}}(h, \bar{u})$ and the flux $\overline{\mathcal{G}}^h = \overline{\mathcal{G}}(h\bar{u}, \phi + gh)$ and $\overline{\mathcal{G}^u} = \overline{\mathcal{G}}(h\bar{u}, \overline{\mathcal{K}}(1, \bar{u}))$ with

$$ (5) \qquad \overline{\mathcal{K}}(h, \psi) = \frac{h}{2} \|\psi\|^2 \quad \text{and} \quad \overline{\mathcal{G}}(\mathcal{Q}, \psi) = \psi \mathcal{Q} $$

More precisely, in case of non-smooth solution, the admissible solutions are defined such that the mechanical energy is decreasing, i.e.

$$ \partial_t (\mathcal{E} + \overline{\mathcal{K}^u}) + \nabla \cdot (\overline{\mathcal{G}}^h + \overline{\mathcal{G}^u}) \le h \partial_t \phi. $$

It is the counterpart of Proposition 1.ii).

### 3.1.2 Entropy-satisfying numerical scheme of (SW)

In the past, several numerical strategies have been proposed that preserve main physical properties, in particular Proposition 2. In the current work, we briefly recall one of them based on an explicit Finite Volume scheme. The scheme of (SW) can be formulate as

$$ (SW^\delta) \qquad U_k^{\text{SW},n+1} = U_k^{\text{SW},n} - \frac{\delta_t^n}{|k|} \sum_{f \in \mathcal{F}_k} \left( \mathcal{F}_f^{\text{SW},n} \cdot \mathcal{N}_k^{k_f} + \mathcal{S}_{f,k}^{\text{SW},n} \right) |f| $$
with $U_k^{\text{sw},n} = (h_k^n, h_k^n \bar{u}_k^n)^t$. Several strategies are described in the literature to compute the numerical flux $\mathcal{F}_f^{\text{sw},n} = (\mathcal{F}_f^{1,n}, \mathcal{F}_f^{2,n})^t$ and the source term $\mathcal{S}_{f,k}^{\text{sw},n}$ of the shallow water model, see [5, 20, 28, 38].

In the following, we ask $(\text{SW}^\delta)$ to satisfy the following properties under a CFL condition like

$$ (6) \qquad \lambda \left( U_k^{\text{sw},n}; U_{\underline{k}_f}^{\text{sw},n} \right) \delta_t^n \le \min \left( \delta_k, \delta_{\underline{k}_f} \right) $$

with $\delta_k$ the compactness of $\underline{k}_f$, i.e. $\delta_k = \frac{|k|}{2}$ if $d=1$ and $\delta_k = \frac{|k|}{\sum_{f \in \mathcal{F}_k} |f|}$ if $d=2$. $\lambda(U_L^{\text{sw}}, U_R^{\text{sw}})$ is an upper limit of the wave speed, depending on the numerical flux $\mathcal{F}^{\text{sw}}$ used, see [5].

**Hypothesis 1.** Under the CFL condition (6) and for any positive initial data, i.e. $h_k^0 \ge 0$, the solution of $(\text{SW}^\delta)$ satisfies

1.i) the water depth positivity, i.e. $h_k^n \ge 0$.

It is the discrete counterpart of Proposition 2.i).

1.ii) the discrete mechanical energy dissipation. More precisely, there exist two flux $\overline{\mathcal{G}}_\delta^h (U_L^{\text{sw}}, U_R^{\text{sw}})$ and $\overline{\mathcal{G}}_\delta^{\overline{u},n} (U_L^{\text{sw}}, U_R^{\text{sw}})$ respectively consistent with $\overline{\mathcal{G}}^h$ and $\overline{\mathcal{G}}^{\overline{u},n}$ given in (5) such that the discrete mechanical energy satisfies

$$ \partial_t^{n+1} (\mathcal{E}_k + \overline{\mathcal{K}}_k^{\overline{u}}) + \frac{1}{|k|} \sum_{f \in \mathcal{F}_k} (\overline{\mathcal{G}}_f^{h,n} + \overline{\mathcal{G}}_f^{\overline{u},n}) \cdot \mathcal{N}_k^{k_f} |f| \le h_k^n \partial_t^{n+1} \phi_k $$

with $\mathcal{E}_k^n = \mathcal{E}(h_k^n)$, $\overline{\mathcal{K}}_k^{\overline{u},n} = \overline{\mathcal{K}}(h_k^n, \overline{u}_k^n)$, $\overline{\mathcal{G}}_f^{h,n} = \overline{\mathcal{G}}_\delta^h (U_k^{\text{sw},n}, U_{\underline{k}_f}^{\text{sw},n})$ and $\overline{\mathcal{G}}_f^{\overline{u},n} = \overline{\mathcal{G}}_\delta^{\overline{u}} (U_k^{\text{sw},n}, U_{\underline{k}_f}^{\text{sw},n})$.

It is the discrete counterpart of Proposition 2.ii).

For homogeneous potential, i.e. $\nabla\phi = 0$, some numerical strategies ensure the entropy stability Hypothesis 1.ii). This is the case, for example, of Godunov solver [21], of the kinetic solver [33] and of the Suliciu relaxation solver [37]. In the case of a source term, the numerical strategy only ensure a semi-discrete entropy stability, i.e. before discretization in time, see [5]. Some implicit or ImEx schemes ensure the entropy stability with the source term, see instead [31], and can also be used with the numerical strategy of the dispersive models presented in the following.

Note that the Hypothesis 1 is not an exhaustive list of properties that we can ask for $(\text{SW})$, for instead the conservation of the steady state. Even if the translation of these properties at the discrete level is an interesting and stimulating question, this is not the purpose of the current work.

## 3.2 Non-hydrostatic model (NH)

### 3.2.1 Governing equations of (NH) and main properties

Let us now focus on the second assumption $H_{yp}^p)$. The simplest model in the hierarchy presented in §2.2 is the non hydrostatic model (NH) which can be formulated
as follows

$$ (NH.a) \qquad \left\{ \begin{array}{l} \partial_t h + \nabla \cdot (h \bar{u}) = 0 \\ \partial_t (h \bar{u}) + \nabla \cdot (h \bar{u} \otimes \bar{u} + \frac{g}{2} h^2 \mathbf{I}_d) + h \nabla \phi = -\nabla \left(h \frac{q_B}{2}\right) - q_B \nabla B \\ \partial_t (h \bar{w}) + \nabla \cdot (h \bar{w} \bar{u}) = q_B \end{array} \right. $$

with the hydrodynamic pressure $q_B$ related to the constraint

$$ (NH.b) \qquad \bar{w} = \partial_t \left(B + \frac{h}{2}\right) + \bar{u} \cdot \nabla \left(B + \frac{h}{2}\right). $$

The initial condition $h(0, x) = h^0(x)$, $\bar{u}(0, x) = \bar{u}^0(x)$ and $\bar{w}(0, x) = \bar{w}^0(x)$ must satisfy the compatibility condition

$$ (8) \qquad \bar{w}^0 = \partial_t B_{|t=0} + \bar{u}^0 \cdot \nabla B - \frac{h^0}{2} \nabla \cdot \bar{u}^0 $$

which corresponds to the vertical-averaged of the compatibility condition (1). For more detail about (NH), we refer to [6]. Let us first recall the main physical properties of (NH).

**Proposition 3.** *Assuming that $h^0(x) \ge 0$, then the solution of (NH) satisfies*

3.i) the water depth positivity, i.e. $h(t, x) \ge 0$.
It is the counterpart of Proposition 1.i).

3.ii) the mechanical energy conservation for a sufficiently smooth solution, i.e.

$$ \partial_t (\mathcal{E} + \bar{\mathcal{K}} \bar{u} + \bar{\mathcal{K}} \bar{w}) + \nabla \cdot (\bar{\mathcal{G}}^h + \bar{\mathcal{G}} \bar{u} + \bar{\mathcal{G}} \bar{w} + \bar{\mathcal{G}}^q)^{\text{NH}} = h \partial_t \phi + q_B \partial_t B $$

with the potential energy $\mathcal{E}(h)$ defined by (2), the kinetic energies $\bar{\mathcal{K}}^\psi = \bar{\mathcal{K}}(h, \psi)$, and the flux $\bar{\mathcal{G}}^h = \bar{\mathcal{G}}(h\bar{u}, \phi + gh)$, $\bar{\mathcal{G}} \bar{u} = \bar{\mathcal{G}}(h\bar{u}, \bar{\mathcal{K}}(1, \bar{u}))$, $\bar{\mathcal{G}} \bar{w} = \bar{\mathcal{G}}(h\bar{u}, \bar{\mathcal{K}}(1, \bar{w}))$, and $\bar{\mathcal{G}}^q = \bar{\mathcal{G}}(h\bar{u}, \frac{q_B}{2})$ defined by (5).
It is the counterpart of Proposition 1.ii).

### 3.2.2 Entropy-satisfying numerical scheme of (NH)

Now we focus on describing a numerical strategy for solving (NH). As in previous works [1, 4, 9, 26], the numerical scheme is based on a splitting strategy between the advection step and the dispersion step. The splitting strategy is a well unknown solution to design numerical scheme for models where a part of the equation is hyperbolic [24]. The advantage of this strategy is to use the classic scheme ($SW^\delta$) for the advection part. In particular, thanks to the formulation (NH), the advection (left-hand side) of the vertical velocity is solved by conventional passive transport scheme. This strategy has already been used in [1] and entropy-satisfying schemes to solve it can be used [5]. However, in the cited works, the dispersion step is solved
using a prediction-correction strategy based on an implicit system on the hydrodynamical pressure $q_B$. Another strategy to solve the dispersive model is to write the equation on horizontal velocity, see [4, 9, 26]. The advantage of this strategy is that dissipative forces such as friction or viscous terms can also be taken into account simultaneously with dispersion. Unfortunately, without the introduction of vertical velocity, this strategy reveals third order derivatives that are not easy to discretize. The scheme proposed in the current work attempts to address the benefits of both strategies. More precisely, the vertical velocity is introduced and advected during the hyperbolic step and the dispersion step is resolved using an elliptical equation on the horizontal velocity. Even if dissipative forces, friction and viscous terms, can be taken into account simultaneously with dispersion as in previous strategies based on horizontal velocity, we do not deal with these terms in this work for readability. We reformulated (NH) as

$$ (9) \qquad \partial_t U^{\text{NH}} + \nabla \cdot F^{\text{NH}}(U^{\text{NH}}) = S^{\text{NH}}(U^{\text{NH}}) + D^{\text{NH}} \\ \text{with} \quad U^{\text{NH}} = \begin{pmatrix} U^{\text{SW}} \\ \hbar \bar{w} \end{pmatrix}, \quad F^{\text{NH}}(U^{\text{NH}}) = \begin{pmatrix} F^{\text{SW}}(U^{\text{SW}}) \\ \hbar \bar{w} \bar{u} \end{pmatrix} \\ S^{\text{NH}}(U^{\text{NH}}) = \begin{pmatrix} S^{\text{SW}}(U^{\text{SW}}) \\ 0 \end{pmatrix} \quad \text{and} \quad D^{\text{NH}} = \begin{pmatrix} 0 \\ -\nabla(h \frac{q_B}{2}) - q_B \nabla B \\ q_B \end{pmatrix}. $$

Let's now assume that the water depth, the horizontal velocity and vertical velocity at the time of iteration $n$ are known, i.e. $h_k^n, \bar{u}_k^n$ and $\bar{w}_k^n$, as for example at the initial condition. Then, the time step can be split into three steps, briefly described as follows:

1. The advection step (NH$^\delta$.a), where the advection and the conservative forces are solved using (SW$^\delta$).

2. The dispersion step (NH$^\delta$.b), where the dispersive operator, and eventually the dissipative forces, are solved using an implicit finite volume method.

### Advection and conservative forces:

This step is based on (SW$^\delta$). More precisely, by neglecting the dispersal operator $D^{\text{NH}}$, the first two unknowns of $U^{\text{NH}}$ are estimated using (SW$^\delta$) while the third is simply advected to the flow as a passive pollutant, see [5]. This strategy was already used for dispersive model in [1]. We write

$$ (NH^\delta \cdot \text{a}) \qquad U_k^{\text{NH},n*} = U_k^{\text{NH},n} - \frac{\delta_t^n}{|k|} \sum_{f \in \mathcal{F}_k} \left( \mathcal{F}_f^{\text{NH},n} \cdot \mathcal{N}_k^{k_f} + \mathcal{S}_{f,k}^{\text{NH},n} \right) |f| $$

with the state vector $U_k^{\text{NH},n} = (U_k^{\text{SW},n}, h_k^n \bar{w}_k^n)^t$, the source term $\mathcal{S}_{f,k}^{\text{NH},n} = (\mathcal{S}_{f,k}^{\text{SW},n}, 0)^t$, the numerical flux

$$ \mathcal{F}_f^{\text{NH},n} \cdot \mathcal{N}_k^{k_f} = \left( \mathcal{F}_f^{\text{SW},n} \cdot \mathcal{N}_k^{k_f} \mathcal{F}_\delta^{\text{up}} \left( \mathcal{F}_f^{1,n} \cdot \mathcal{N}_k^{k_f} \left( \bar{w}_k^n, \bar{w}_{k_f}^n \right) \right) \right) $$
with the up-wind scheme

$$ (11) \qquad \mathcal{F}_{\delta}^{\mathrm{up}}(Q, (\psi_L, \psi_R)) = \psi_L Q_+ - \psi_R Q_- $$

and the positive and negative part functions are defined by $\psi_{\pm} = \frac{|\psi| \pm \psi}{2}$.

### Dispersion and dissipative forces:

Now let us introduce the scheme of the dispersion step. Since the water depth is not affected by this step, we have defined $h_k^{n+1} = h_k^{n*}$. Then, wherever the water depth does not vanishes, the horizontal velocity is corrected by the following implicit linear scheme

$$ (NH^{\delta}.b) \qquad \alpha_k^{\mathrm{NH}} \bar{u}_k^{n+1} + \nabla_k^{\delta} (\mu^{\mathrm{NH}} \cdot \bar{u}^{n+1}) - \mu_k^{\mathrm{NH}} \nabla_k^{\delta} \bar{u}^{n+1} - \nabla_k^{\delta} (\kappa^{\mathrm{NH}} \nabla_k^{\delta} \bar{u}^{n+1}) = \beta_k^{\mathrm{NH}} $$

with the discrete operators $\nabla_k^{\delta}$ and $\nabla_k^{\delta}$, defined by (3) and the parameters

$$ \begin{gathered} \alpha_k^{\mathrm{NH}} = h_k^{n*} (\mathbf{I}_d + \nabla_k^{\delta} B^{n+1} \otimes \nabla_k^{\delta} B^{n+1}), \\ \mu_k^{\mathrm{NH}} = \frac{(h_k^{n*})^2}{2} \nabla_k^{\delta} B^{n+1}, \qquad \kappa_k^{\mathrm{NH}} = \frac{(h_k^{n*})^3}{4}, \\ \text{and} \qquad \beta_k^{\mathrm{NH}} = h_k^{n*} \bar{u}_k^{n*} + h_k^{n*} (\bar{w}_k^{n*} - \partial_t^{n+1} B_k) \nabla_k^{\delta} B^{n+1} \\ \qquad + \nabla_k^{\delta} \left( \frac{(h_k^{n*})^2}{2} (\bar{w}^{n*} - \partial_t^{n+1} B) \right). \end{gathered} $$

In the dry zone, the scheme leads to a trivial equation and the resulting linear system is not well-posed. To overcome this issue, in all dry cells, $(NH^{\delta}.b)$ is replaced by $\bar{u}_k^{n+1} = 0$. The solution is obviously influenced by the choice of velocity in dry areas. However, this choice is motivated by the conservation of the mechanical energy at the dry front, see Proposition 4. More precisely, thanks to this choice, the energy flux of the dispersive part vanishes on the dry front.

Finally, the vertical velocity is reconstructed using the following scheme

$$ (NH^{\delta}.c) \qquad \bar{w}_k^{n+1} = \partial_t^{n+1} B_k + \bar{u}_k^{n+1} \cdot \nabla_k^{\delta} B^{n+1} - \frac{h_k^{n+1}}{2} \nabla_k^{\delta} \bar{u}^{n+1} $$

where the discrete time derivative is defined by (4) and the discrete operators are defined by (3).

Let's explain how the scheme is obtained. $(NH^{\delta}.c)$ is obtained after a rewriting of $(NH.b)$ using the mass conservation (first equation of (NH.a)) as

$$ \bar{w} = \partial_t B + \bar{u} \cdot \nabla B - \frac{h}{2} \nabla \cdot \bar{u}. $$

Even if other writings of the constraint can be considered, a discretization of this one ensures the entropy stability at the discret level, see Proposition 4. $(NH^{\delta}.b)$ is obtained by considering the following discretization of the dispersive operator

$$ (13) \qquad U_k^{\mathrm{NH},n+1} = U_k^{\mathrm{NH},n*} + \delta_t^n \mathcal{D}_k^{\mathrm{NH},n+1} (q_B^{n+1}) $$
with

$$ \mathcal{D}_k^{\text{NH},n}(q_B) = \begin{pmatrix} 0 \\ -\nabla_k^\delta \left(h^n \frac{q_B}{2}\right) - q_{B,k} \nabla_k^\delta B^n \end{pmatrix}. $$

Together with (NH$^\delta$.c), they form a $(d+2) \times N$ linear system (without $h$). Fortunately, this system can be reduced by simple substitution and leads to the $d \times N$ linear system (NH$^\delta$.b).

**Proposition 4.** *Assuming Hypothesis 1, $h_k^0 \ge 0$ and the CFL condition (6) holds. Then the solution of (NH$^\delta$) satisfies*

4.i) the water depth positivity, i.e. $h_k^n \ge 0$.

It is the discrete counterpart of Proposition 3.i).

4.ii) the discrete mechanical energy dissipation, i.e.

$$ \partial_t^{n+1} (\mathcal{E}_k + \overline{\mathcal{K}}_k^{\overline{u}} + \overline{\mathcal{K}}_k^{\overline{w}}) + \frac{1}{|k|} \sum_{f \in \mathcal{F}_k} \left( \overline{\mathcal{G}}_f^{h,n} + \overline{\mathcal{G}}_f^{\overline{u},n} + \overline{\mathcal{G}}_f^{\overline{w},n} + \overline{\mathcal{G}}_f^{q^{\text{NH},n+1}} \right) \cdot \mathcal{N}_k^{k_f} |f| \le h_k^n \partial_t^{n+1} \phi_k + q_{B,k}^{n+1} \partial_t^{n+1} B_k $$

with the energies $\mathcal{E}_k^n = \mathcal{E}(h_k^n)$, $\overline{\mathcal{K}}_k^{\psi,n} = \overline{\mathcal{K}}(h_k^n, \psi_k^n)$ defined by (5),

the flux $\overline{\mathcal{G}}_f^{h,n} = \overline{\mathcal{G}}_\delta^h \left( \left( h_k^n, h_k^n \overline{u}_k^n \right), \left( h_{\overline{k}_f}^n, h_{\overline{k}_f}^n \overline{u}_{\overline{k}_f}^n \right) \right)$, $\overline{\mathcal{G}}_f^{\overline{u},n} = \overline{\mathcal{G}}_\delta^{\overline{u}} \left( \left( h_k^n, h_k^n \overline{u}_k^n \right), \left( h_{\overline{k}_f}^n, h_{\overline{k}_f}^n \overline{u}_{\overline{k}_f}^n \right) \right)$ introduced in Hypothesis 1.ii),

$\overline{\mathcal{G}}_f^{\overline{w},n} = \mathcal{F}_\delta^{up} \left( \mathcal{F}_f^{1,n} \cdot \mathcal{N}_k^{k_f}, \left( \frac{(\overline{w}_k^n)^2}{2}, \frac{(\overline{w}_{\overline{k}_f}^n)^2}{2} \right) \right)$ defined

by (11) and $\overline{\mathcal{G}}_f^{q^{\text{NH},n}} = \overline{\mathcal{G}}_\delta^{\overline{q}} \left( \left( \frac{h_k^n}{\overline{u}_k^n}, \frac{h_{\overline{k}_f}^n}{\overline{u}_{\overline{k}_f}^n} \right), \left( \frac{q_{B,k}^n}{2}, \frac{q_{B,k_f}^n}{2} \right) \right)$ defined by

$$ \overline{\mathcal{G}}_\delta^{\overline{q}} \left( \left( \frac{h_L}{\overline{q}_L}, \frac{h_R}{\overline{q}_R} \right), \left( \frac{u_L}{\overline{q}_L}, \frac{u_R}{\overline{q}_R} \right) \right) = \frac{h_L \overline{q}_L + h_R \overline{q}_R}{2} \frac{\overline{u}_L + \overline{u}_R}{2} - \frac{h_R \overline{q}_R - h_L \overline{q}_L}{2} \frac{\overline{u}_R - \overline{u}_L}{2}. $$

It is the discrete counterpart of Proposition 3.ii).

*Proof.* Proposition 4.i) is an obvious consequence of Hypothesis 1.i).

Let's focus on mechanical energy conservation Proposition 4.ii). Thanks to Hypothesis 1.ii), we have

$$ \begin{align} \mathcal{E}_k^{n\star} + \overline{\mathcal{K}}_k^{\overline{u},n\star} + \frac{\delta_t^n}{|k|} \sum_{f \in \mathcal{F}_k} \left( \overline{\mathcal{G}}_f^{h,n} + \overline{\mathcal{G}}_f^{\overline{u},n} \right) \cdot \mathcal{N}_k^{k_f} |f| \nonumber \\ \le \mathcal{E}_k^n + \overline{\mathcal{K}}_k^{\overline{u},n} + h_k^n (\phi_k^{n+1} - \phi_k^n) \end{align} $$
with $\mathcal{E}_k^{n\star} = \mathcal{E}(h_k^{n\star})$ and $\overline{\mathcal{K}}_k^{\psi,n\star} = \overline{\mathcal{K}}(h_k^{n\star}, \psi_k^{n\star})$.

Similarly, a classical result of the Suliciu relaxation is the following entropy inequality as long as the positivity is ensured, see [5, Chapter 2.7]

$$ (15) \qquad \overline{\mathcal{K}}_k^{\overline{w},n\star} + \frac{\delta_t^n}{|k|} \sum_{f \in \mathbb{F}_k} \overline{\mathcal{G}}_f^{\overline{w},n} |f| \le \overline{\mathcal{K}}_k^{\overline{w},n}. $$

Now, we focus on the dispersion step (NH$^\delta$.b)-(NH$^\delta$.c). As explained above, the dispersion step (NH$^\delta$.b) is equivalent to (13). By multiplying the second equation of (13) by $\overline{u}_k^{n+1}$, it yields

$$ (16) \qquad \overline{\mathcal{K}}_k^{\overline{u},n+1} \le \overline{\mathcal{K}}_k^{\overline{u},n\star} - \delta_t^n \left( \overline{u}_k^{n+1} \cdot \nabla_k^\delta \left( h^{n+1} \frac{q_B^{n+1}}{2} \right) + q_{B,k}^{n+1} \overline{u}_k^{n+1} \cdot \nabla_k^\delta B^{n+1} \right). $$

Thanks to the centered operator (3), the following relation holds

$$ \overline{u}_k^{n+1} \cdot \nabla_k^\delta \left( h^{n+1} \frac{q_B^{n+1}}{2} \right) = \frac{1}{|k|} \sum_{f \in \mathbb{F}_k} \overline{\mathcal{G}}_f^{q^{NH},n+1} \cdot \mathcal{N}_k^{\frac{k}{f}} |f| - q_{B,k}^{n+1} \frac{h_k^{n+1}}{2} \nabla_k^\delta \cdot \overline{u}^{n+1}. $$

Then by multiplying the third equation of (13) by $\overline{w}_k^{n+1}$, it yields

$$ (17) \qquad \overline{\mathcal{K}}_k^{\overline{w},n+1} \le \overline{\mathcal{K}}_k^{\overline{w},n\star} + \delta_t^n q_{B,k}^{n+1} \overline{w}_k^{n+1}. $$

Summing the energy dissipations (16) and (17), it yields

$$ (18) \qquad \overline{\mathcal{K}}_k^{\overline{u},n+1} + \overline{\mathcal{K}}_k^{\overline{w},n+1} + \frac{\delta_t^n}{|k|} \sum_{f \in \mathbb{F}_k} \overline{\mathcal{G}}_f^{q^{NH},n+1} \cdot \mathcal{N}_k^{\frac{k}{f}} |f| \le \overline{\mathcal{K}}_k^{\overline{u},n\star} + \overline{\mathcal{K}}_k^{\overline{w},n\star} \\ + \delta_t^n q_{B,k}^{n+1} \left( \overline{w}_k^{n+1} - \overline{u}_k^{n+1} \cdot \nabla_k^\delta B^{n+1} + \frac{h_k^{n+1}}{2} \nabla_k^\delta \cdot \overline{u}^{n+1} \right). $$

Finally, we conclude by using (NH$^\delta$.c). At wet/dry fronts, the water level vanishes in a cell, for instance $h_R = 0$. Then the energy flux is proportional to the velocity in the dry cell, i.e. $\overline{\mathcal{G}}_\delta^q = \frac{h_L \overline{q}_L}{2} \overline{u}_R$. By setting $\overline{u}_R = 0$, the flux through the wet/dry fronts vanishes. $\square$

**Remark 1.** The dissipative term such are friction and viscosity can easily be treated simultaneously with dispersion. More precisely, it only adds the corresponding reaction and diffusion terms in (NH$^\delta$.b).

**Remark 2.** A drawback of (NH$^\delta$.b) is its large stencil, i.e. the 5-point stencil in 1D. Some may want to replace (NH$^\delta$.b) with a more compact and simpler discrete operator, see Appendix A. However, the entropy stability Proposition 4.ii) is lost and in practice the numerical results are less accurate, see 4.2.
**Remark 3.** Some people may rather like to solve the linear system by keeping one of the other unknowns after substitution. For example, in [1], the dispersion step is solved by maintaining the hydrodynamic pressure $q_B$ as unknown to the linear system. This choice is particularly attractive for multidimensional simulations, since $\bar{u} \in \mathbb{R}^d$ whereas $q_B \in \mathbb{R}$. However, this formulation is less practical for imposing zero velocity in dry areas. Moreover, it is not possible to treat friction and viscosity, if any, in the same step, contrary to the formulation (NH$^\delta$.b).

Let us define the set of admissible solutions for (NH)

$$ (19) \quad \mathcal{A}_r^{\text{NH},n} = \left\{ (\bar{u}, \bar{w}) \in (\mathbb{R}^N)^{d+1} \mid \left\| \partial_t^n B + \bar{u}_k \cdot \nabla^\delta B^n - \frac{h^n}{2} \nabla^\delta \cdot \bar{u} - \bar{w} \right\|_{\delta} \le r \right\}. $$

with $\|\cdot\|_\delta$ is a discret norm on the tessellation $\mathbb{T}$. In the following we use the $L^2$-norm

$$ \|\phi\|_{\delta} = \sum_{k \in \mathbb{T}} |k| \phi_k^2. $$

**Proposition 5.** The dispersion step (NH$^\delta$.b) is a projection on the set of velocity satisfying the constrain (NH$^\delta$.c).

More precisely $(\bar{u}^{n+1}, \bar{w}^{n+1}) \in \mathcal{A}_0^{\text{NH},n+1}$

and if $(\bar{u}^{n*}, \bar{w}^{n*}) \in \mathcal{A}_0^{\text{NH},n+1}$ then $(\bar{u}^{n+1}, \bar{w}^{n+1}) = (\bar{u}^{n*}, \bar{w}^{n*})$.

## 3.3 Serre-Green-Naghdi model (GN)

### 3.3.1 Governing equations of (GN) and main properties

As it was mentioned in [17], the non-hydrostatic model (NH) is not completely coherent with (E). More precisely, the vectorial space of the approximation for $u$ and $w$, i.e. both in $\mathbb{P}_0(z)$, is not coherent with the divergence free condition. To make the approximation coherent with the divergence free condition when $u$ is in $\mathbb{P}_0(z)$, $w$ must be approximated in $\mathbb{P}_1(z)$, see §2.2.3 and [17] for more details. It is the case for the Serre-Green-Naghdi model (GN), see [22, 32, 35].

To propose a scheme in the spirit of (NH$^\delta$), (GN) is reformulated as follows

$$ (\text{GN}.a) \quad \left\{ \begin{array}{l} \partial_t h + \nabla \cdot (h \bar{u}) = 0 \\ \partial_t (h \bar{u}) + \nabla \cdot (h \bar{u} \otimes \bar{u} + \frac{g}{2} h^2 \mathbf{I}_d) = -h \nabla \phi - \nabla (h \bar{q}) - q_B \nabla B \\ \partial_t (h \bar{w}) + \nabla \cdot (h \bar{w} \bar{u}) = q_B \\ \partial_t (h \sigma) + \nabla \cdot (h \sigma \bar{u}) = \sqrt{3} (2 \bar{q} - q_B) \end{array} \right. $$

with the hydrodynamic pressures $\bar{q}$ and $q_B$ related to constraints (NH.b) and

$$ (\text{GN}.b) \qquad \sigma = \frac{\partial_t h + \bar{u} \cdot \nabla h}{2\sqrt{3}}. $$

The initial condition $h(0, x) = h^0(x)$, $\bar{u}(0, x) = \bar{u}^0(x)$, $\bar{w}(0, x) = \bar{w}^0(x)$ and $\sigma(0, x) = \sigma^0(x)$ must satisfy the compatibility conditions (8) and

$$ (21) \qquad \sigma^0 = -\frac{h^0}{2\sqrt{3}} \nabla \cdot \bar{u}^0 $$
which corresponds to the vertical-averaged of the standard deviation of the compatibility condition (1). In [17] it is shown that (GN) is equivalent, at least for smooth solution, to the classical Peregrine formulation [32]. Let us first recall the main physical properties of (GN).

**Proposition 6.** *Assuming that $h^0(x) \ge 0$, then the solution of (GN) satisfies*

6.i) the water depth positivity, i.e. $h(t, x) \ge 0$.
It is the counterpart of Proposition 1.i).

6.ii) the mechanical energy conservation for a sufficiently smooth solution, i.e.

$$ \partial_t (\mathcal{E} + \overline{\mathcal{K}}^u + \overline{\mathcal{K}}^w + \overline{\mathcal{K}}^\sigma) + \nabla \cdot (\overline{\mathcal{G}}^h + \overline{\mathcal{G}}^u + \overline{\mathcal{G}}^w + \overline{\mathcal{G}}^\sigma + \overline{\mathcal{G}}^q) = h \partial_t \phi + q_B \partial_t B $$

with the potential energy $\mathcal{E}$ (h) defined by (2), the kinetic energies $K^\psi = \overline{\mathcal{K}}(h, \psi)$ and the flux $\overline{\mathcal{G}}^h = \overline{\mathcal{G}}(h\bar{u}, \phi + gh)$, $\overline{\mathcal{G}}^u = \overline{\mathcal{G}}(h\bar{u}, \overline{\mathcal{K}}(1, \bar{u}))$, $\overline{\mathcal{G}}^w = \overline{\mathcal{G}}(h\bar{u}, \overline{\mathcal{K}}(1, \bar{w}))$, $\overline{\mathcal{G}}^\sigma = \overline{\mathcal{G}}(h\bar{u}, \overline{\mathcal{K}}(1, \sigma))$, $\overline{\mathcal{G}}^q = \overline{\mathcal{G}}(h\bar{u}, \bar{q})$ defined by (5).
It is the counterpart of Proposition 1.ii).

### 3.3.2 Entropy-satisfying numerical scheme of (GN)

Now we focus on describing a numerical strategy for solving (GN). Similarly to (NH$^\delta$), (GN.a) is written as

$$ \begin{gathered} \partial_t U^{\text{GN}} + \nabla \cdot F^{\text{GN}}(U^{\text{GN}}) = S^{\text{GN}}(U^{\text{GN}}) + D^{\text{GN}} \\ \text{with} \quad U^{\text{GN}} = \begin{pmatrix} U^{\text{NH}} \\ h\sigma \end{pmatrix}, \quad F^{\text{GN}}(U^{\text{GN}}) = \begin{pmatrix} F^{\text{NH}}(U^{\text{NH}}) \\ h\sigma u \end{pmatrix} \\ S^{\text{GN}}(U^{\text{GN}}) = \begin{pmatrix} S^{\text{NH}}(U^{\text{NH}}) \\ 0 \end{pmatrix} \quad \text{and} \quad D^{\text{GN}} = \begin{pmatrix} 0 \\ -\nabla(h\bar{q}) - q_B \nabla B \\ q_B \\ \sqrt{3}(2\bar{q} - q_B) \end{pmatrix} \end{gathered} $$

and the numerical strategy is based on the same two steps: advection and conservative forces, dispersion and dissipative forces.

**Advection and conservative forces:**

As done for (NH$^\delta$.a), the advection and conservative forces are realized using the strategy (SW$^\delta$) with an up-wind scheme for the vertical velocities $\bar{w}_k^n$ and $\sigma_k^n$. We write

$$ (GN^\delta.a) \qquad U_k^{\text{GN},n*} = U_k^{\text{GN},n} - \frac{\delta_t^n}{|k|} \sum_{f \in \mathcal{F}_k} \left( \mathcal{F}_f^{\text{GN},n} \cdot \mathcal{N}_k^{k_f} + \mathcal{S}_{f,k}^{\text{GN},n} \right) |f| $$
with the state vectors $U_k^{\text{GN},n} = (U_k^{\text{NH},n}, h_k^n \sigma_k^n)^t$, the numerical flux

$$ \mathcal{F}_f^{\text{GN},n} \cdot \mathcal{N}_k^{\frac{k}{f}} = \left( \begin{array}{c} \mathcal{F}_f^{\text{NH},n} \cdot \mathcal{N}_k^{\frac{k}{f}} \\ \mathcal{F}_\delta^{\text{up}} (\mathcal{F}_f^{1,n} \cdot \mathcal{N}_k^{\frac{k}{f}}, (\sigma_k^n, \sigma_{\frac{k}{f}}^n)) \end{array} \right) $$

with the upwind flux defined by (11) and the source term $\mathcal{S}_{f,k}^{\text{GN},n} = (\mathcal{S}_{f,k}^{\text{NH},n}, 0)^t$.

**Dispersion and dissipative forces:**

As for $(NH^\delta)$, we set $h_k^{n+1} = h_k^{n*}$. In wet cell, i.e. $h_k^{n+1} > 0$, we write

$$ (GN^\delta.\text{b}) \qquad \alpha_k^{\text{GN}} \bar{u}_k^{n+1} + \nabla_k^\delta (\mu_k^{\text{GN}} \cdot \bar{u}^{n+1}) - \mu_k^{\text{GN}} \nabla_k^\delta \bar{u}^{n+1} - \nabla_k^\delta (\kappa_k^{\text{GN}} \nabla_k^\delta \cdot \bar{u}^{n+1}) = \beta_k^{\text{GN}} $$

where the parameters read

$$ \alpha_k^{\text{GN}} = \alpha_k^{\text{NH}}, \quad \mu_k^{\text{GN}} = \mu_k^{\text{NH}}, \quad \kappa_k^{\text{GN}} = \frac{4}{3} \kappa_k^{\text{NH}} \quad \text{and} \quad \beta_k^{\text{GN}} = \beta_k^{\text{NH}} + \frac{1}{2\sqrt{3}} \nabla_k^\delta ((h^{n*})^2 \sigma^{n*}). $$

and in dry cell we set $\bar{u}_k^{n+1} = 0$. Finally, the vertical velocity is reconstructed using $(NH^\delta.\text{c})$ and

$$ (GN^\delta.\text{c}) \qquad \sigma_k^{n+1} = -\frac{h_k^{n+1}}{2\sqrt{3}} \nabla_k^\delta \bar{u}^{n+1}. $$

The scheme $(GN^\delta.\text{b})$ was obtained similarly to $(NH^\delta.\text{b})$. More precisely, the dispersive operator is discretized by

$$ (23) \qquad U_k^{\text{GN},n+1} = U_k^{\text{GN},n*} + \delta_t^n \mathcal{D}_k^{\text{GN},n+1} (\bar{q}^{n+1}, \bar{q}_B^{n+1}) $$

with

$$ \mathcal{D}_k^{\text{GN},n} (\bar{q}, \bar{q}_B) = \begin{pmatrix} 0 \\ -\nabla_k^\delta (h^n \bar{q}) - q_{B,k} \nabla_k^\delta B^n \\ q_{B,k} \\ \sqrt{3} (2 \bar{q}_k - q_{B,k}) \end{pmatrix}. $$

Together with $(GN^\delta.\text{c})$, it forms a $(d+3) \times N$ linear system which can be reduced by simple substitution and leads to the $d \times N$ linear system $(GN^\delta.\text{b})$.

**Proposition 7.** *Assuming Hypothesis 1, $h_k^0 \ge 0$ and the CFL condition (6) holds. Then the solution of $(GN^\delta)$ satisfies*

7.i) the water depth positivity, i.e. $h_k^n \ge 0$. It is the discrete counterpart of Proposition 6.i).
7.ii) the discrete mechanical energy dissipation, i.e.

$$
\begin{align*}
& \partial_t^{n+1} (\mathcal{E}_k + \overline{\mathcal{K}}_k \overline{\bar{u}} + \overline{\mathcal{K}}_k \overline{\bar{w}} + \overline{\mathcal{K}}_k \overline{\sigma}) \\
& \qquad + \frac{1}{|k|} \sum_{f \in \mathcal{F}_k} (\overline{\mathcal{G}}_f^{h,n} + \overline{\mathcal{G}}_f^{\bar{u},n} + \overline{\mathcal{G}}_f^{\bar{w},n} + \overline{\mathcal{G}}_f^{\sigma,n} + \overline{\mathcal{G}}_f^{q^{\text{GN},n+1}}) \cdot \mathcal{N}_k^{\underline{k}_f} |f| \\
& \leq h_k^n \partial_t^{n+1} \phi_k + q_{B,k}^{n+1} \partial_t^{n+1} B_k
\end{align*}
$$

with the energies $\mathcal{E}_k^n = \mathcal{E}(h_k^n)$, $\overline{\mathcal{K}}_k^{\psi,n} = \overline{\mathcal{K}}(h_k^n, \psi_k^n)$ defined by (5),
the flux $\overline{\mathcal{G}}_f^{h,n} = \overline{\mathcal{G}}_\delta^h \left( \left( h_k^n \overline{u}_k^n \right), \left( h_{\underline{k}_f}^n \overline{u}_{\underline{k}_f}^n \right) \right)$, $\overline{\mathcal{G}}_f^{\bar{u},n} = \overline{\mathcal{G}}_\delta^{\bar{u}} \left( \left( h_k^n \overline{u}_k^n \right), \left( h_{\underline{k}_f}^n \overline{u}_{\underline{k}_f}^n \right) \right)$ intro-

duced in Hypothesis 1.ii), $\overline{\mathcal{G}}_f^{\bar{w},n} = \mathcal{F}_\delta^{up} \left( \mathcal{F}_f^{1,n} \cdot \mathcal{N}_k^{\underline{k}_f}, \left( \frac{(\overline{w}_k^n)^2}{2}, \frac{(\overline{w}_{\underline{k}_f}^n)^2}{2} \right) \right)$, $\overline{\mathcal{G}}_f^{\sigma,n} =$
$\mathcal{F}_\delta^{up} \left( \mathcal{F}_f^{1,n} \cdot \mathcal{N}_k^{\underline{k}_f}, \left( \frac{(\sigma_k^n)^2}{2}, \frac{(\sigma_{\underline{k}_f}^n)^2}{2} \right) \right)$ defined by (11) and $\overline{\mathcal{G}}_f^{q^{\text{GN},n}} = \overline{\mathcal{G}}_\delta^{\overline{q}} \left( \left( \frac{h_k^n}{\overline{u}_k^n}, \frac{h_{\underline{k}_f}^n}{\overline{u}_{\underline{k}_f}^n} \right), \left( \frac{\overline{q}_k^n}{\overline{q}_k^n}, \frac{\overline{q}_{\underline{k}_f}^n}{\overline{q}_{\underline{k}_f}^n} \right) \right)$

is introduced in Proposition 4.ii).

It is the discrete counterpart of Proposition 6.ii).

Proof. Proposition 7.i) is an obvious consequence of Hypothesis 1.i).

Let us focus on the mechanical energy conservation Proposition 7.ii). The proof is
similar to the one of Proposition 4.ii). Thanks to Hypothesis 1.ii), we obtain (14).
Using the properties of the up-wind scheme, we obtain (15) and

$$
(24) \qquad \overline{\mathcal{K}}_k^{\sigma,n*} + \frac{\delta_t^n}{|k|} \sum_{f \in \mathcal{F}_k} \overline{\mathcal{G}}_f^{\sigma,n} |f| \le \overline{\mathcal{K}}_k^{\sigma,n}.
$$

Now multiplying the second equation of (23) by $\overline{u}_k^{n+1}$, it reads

$$
\begin{equation}
\begin{aligned}
& \overline{\mathcal{K}}_k^{\bar{u},n+1} + \frac{\delta_t^n}{|k|} \sum_{f \in \mathcal{F}_k} \overline{\mathcal{G}}_f^{q^{\text{GN}},n+1} \cdot \mathcal{N}_k^{\underline{k}_f} |f| \le \overline{\mathcal{K}}_k^{\bar{u},n*} \\
& \qquad + \delta_t^n \left( \overline{q}_k^{n+1} h_k^{n+1} \nabla_k^\delta \cdot \bar{u}^{n+1} - q_{B,k}^{n+1} \bar{u}_k^{n+1} \cdot \nabla_k^\delta B^{n+1} \right).
\end{aligned}
\tag{25}
\end{equation}
$$

Then multiplying the third equation of (23) by $\overline{w}_k^{n+1}$, we get (17) and the fourth by $\sigma_k^{n+1}$, we get

$$
(26) \quad \overline{\mathcal{K}}_k^{\sigma,n+1} \le \overline{\mathcal{K}}_k^{\sigma,n*} + \delta_t^n \sqrt{3} (2\overline{q}_k^{n+1} - q_{B,k}^{n+1}) \sigma_k^{n+1}
$$

Summing the energy dissipations (25), (17) and (26), it yields

$$
\begin{align}
& \overline{\mathcal{K}}_k^{\bar{u},n+1} + \overline{\mathcal{K}}_k^{\bar{w},n+1} + \overline{\mathcal{K}}_k^{\sigma,n+1} + \frac{1}{|k|} \sum_{f \in \mathcal{F}_k} \overline{\mathcal{G}}_f^{q^{\text{GN}},n+1} \cdot \mathcal{N}_k^{\underline{k}_f} |f| \nonumber \\
& \le \overline{\mathcal{K}}_k^{\bar{u},n*} + \overline{\mathcal{K}}_k^{\bar{w},n*} + \overline{\mathcal{K}}_k^{\sigma,n*} + \delta_t^n 2\sqrt{3}\overline{q}_k^{n+1} \left( \sigma_k^{n+1} + \frac{h_k^{n+1}}{2\sqrt{3}} \nabla_k^\delta \cdot \overline{u}^{n+1} \right) \nonumber \\
& \quad + \delta_t^n q_{B,k}^{n+1} \left( \overline{w}_k^{n+1} - \left( \overline{u}_k^{n+1} \cdot \nabla_k^\delta B^{n+1} + \sqrt{3}\sigma_k^{n+1} \right) \right). \tag{27}
\end{align}
$$
We conclude using (GN$^{δ}$.c) and summing with (14), (15) and (24).

<figure><img src="image_2.png" /></figure>

**Remark 4.** Remarks 1 and 2 still holds for (GN$^{δ}$). Remark 3 still holds for (GN$^{δ}$) with a slight revision of the conclusion. The system could be solved either with $\bar{u} \in \mathbb{R}^d$ or with $(\bar{q}, q_B) \in \mathbb{R}^2$ as the main unknown. The first solution is clearly preferable while $d=1$ in addition to the fact that the dry areas are treated simply.

Let us define the set of admissible solutions for (GN)

$$ (28) \qquad A_r^{\text{GN},n} = \left\{ (\bar{u}, \bar{w}, \sigma) \in A_r^{\text{NH},n} \times \mathbb{R}^N \;\middle|\; \left\| \sigma + \frac{h^n}{2\sqrt{3}} \nabla^\delta \cdot \bar{u} \right\|_\delta \le r \right\}. $$

**Proposition 8.** The second step of the scheme is actually a projection on the set of velocity satisfying the constrain (NH$^{δ}$.c) and (GN$^{δ}$.c).

More precisely, $(\bar{u}^{n+1}, \bar{w}^{n+1}, \sigma^{n+1}) \in A_0^{\text{GN},n+1}$ and if $(\bar{u}^{n*}, \bar{w}^{n*}, \sigma^{n*}) \in A_0^{\text{GN},n+1}$ then $(\bar{u}^{n+1}, \bar{w}^{n+1}, \sigma^{n+1}) = (\bar{u}^{n*}, \bar{w}^{n*}, \sigma^{n*})$.

## 3.4 Boundary conditions treatment

In this section, a simple treatment of boundary conditions is presented. Even at the continuous level, it is not clear what boundary conditions are compatible with the dispersive models. We only present here the numerical treatment of boundary conditions, i.e. how they are prescribed in practice. It is not clear that the following strategy is fully consistent with the ongoing framework and a rigorous analysis is required.

Let $\partial\mathbb{F}$ be the set of face at the boundary and for any $\tilde{f} \in \partial\mathbb{F}$ let $\tilde{k} \in \mathbb{T}$ be in the computational domain whereas its neighbor $\tilde{k}_{\tilde{f}}$ is the ghost cell.

The advection step consists of the resolution of a system of hyperbolic equations. For a (not exhaustive) list of boundary conditions for (SW), see [20, Chapter V]. The dispersive models also required a boundary condition for $\bar{w}$ and $\sigma$ when the flow is incoming, i.e. $\mathcal{F}^{1,n} \cdot \mathcal{N}_{\tilde{k}}^{\tilde{k}_{\tilde{f}}} < 0$. In practice, we impose a value, usually zero, a the ghost cell.

Let us now focus on the dispersive step (GN$^{δ}$.b) (resp. (NH$^{δ}$.b)). Even if the horizontal velocity $\bar{u}_{\tilde{k}_{\tilde{f}}}^{n+1}$ is already known, several other unknown are required to write the scheme (GN$^{δ}$.b) (resp. (NH$^{δ}$.b)) in $\tilde{k}$. From (23) (resp. (13)), it is clear that only $\bar{q}_{\tilde{k}_{\tilde{f}}}^{n+1}$ (resp. $q_{B,\tilde{k}_{\tilde{f}}}^{n+1}$) is actually required. The strategy to impose the hydrodynamic pressure at the boundary was already propose in [1]. In this paper, we only consider the Neumann condition $\bar{q}_{\tilde{k}_{\tilde{f}}}^{n+1} = \bar{q}_{\tilde{k}}^{n+1}$ (resp. $q_{B,\tilde{k}_{\tilde{f}}}^{n+1} = q_{B,\tilde{k}}^{n+1}$). More precisely, the
scheme in the boundary cells $\tilde{k}$ reads

$$ \alpha_{\tilde{k}}^{\text{XX}} \tilde{u}_{\tilde{k}}^{n+1} + \frac{1}{2|\tilde{k}|} \left( \sum_{f \in \mathbb{F}_{\tilde{k}} \setminus \partial \mathbb{F}} \mu_{\tilde{k}_f}^{\text{XX}} \cdot \overline{u}_{\tilde{k}_f}^{n+1} \mathcal{N}_{\tilde{k}}^{\tilde{k}_f} |f| + \sum_{\tilde{f} \in \mathbb{F}_{\tilde{k}} \cap \partial \mathbb{F}} \mu_{\tilde{k}_f}^{\text{XX}} \cdot \overline{u}_{\tilde{k}_f}^{n+1} \mathcal{N}_{\tilde{k}}^{\tilde{k}_f} |\tilde{f}| \right) - \mu_{\tilde{k}}^{\text{XX}} \nabla_{\tilde{k}}^{\delta} \cdot \overline{u}^{n+1} \\ - \frac{1}{2|\tilde{k}|} \left( \sum_{f \in \mathbb{F}_{\tilde{k}} \setminus \partial \mathbb{F}} \kappa_{\tilde{k}_f}^{\text{XX}} \nabla_{\tilde{k}_f}^{\delta} \cdot \overline{u}^{n+1} \mathcal{N}_{\tilde{k}}^{\tilde{k}_f} |f| + \sum_{\tilde{f} \in \mathbb{F}_{\tilde{k}} \cap \partial \mathbb{F}} \kappa_{\tilde{k}_f}^{\text{XX}} \nabla_{\tilde{k}_f}^{\delta} \cdot \overline{u}^{n+1} \mathcal{N}_{\tilde{k}}^{\tilde{k}_f} |\tilde{f}| \right) = \tilde{\beta}_{\tilde{k}}^{\text{XX}} $$

with $XX = NH$ for $(NH^\delta)$, the parameters defined in $(NH^\delta.b)$ and

$$ \mu_{\tilde{k}}^{\text{NH}} = \frac{h_{\tilde{k} \tilde{f}}^{n+1} h_{\tilde{k}}^{n+1}}{2} \nabla_{\tilde{k}}^{\delta} B^{n+1}, \quad \kappa_{\tilde{k}}^{\text{NH}} = \frac{h_{\tilde{k} \tilde{f}}^{n+1} (h_{\tilde{k}}^{n+1})^2}{4}, $$

$$ \tilde{\beta}_{\tilde{k}}^{\text{NH}} = h_{\tilde{k}}^{n+1} \overline{u}_{\tilde{k}}^{n\star} + h_{\tilde{k}}^{n+1} (\overline{w}_{\tilde{k}}^{n\star} - \partial_t^{n+1} B_{\tilde{k}}) \nabla_{\tilde{k}}^{\delta} B^{n+1} \\ + \frac{1}{4|\tilde{k}|} \left( \sum_{f \in \mathbb{F}_{\tilde{k}} \setminus \partial \mathbb{F}} \left( h_{\tilde{k}_f}^{n+1} \right)^2 \left( \overline{w}_{\tilde{k}_f}^{n\star} - \partial_t^{n+1} B_{\tilde{k}_f} \right) \mathcal{N}_{\tilde{k}}^{\tilde{k}_f} |f| \right. \\ \left. + \sum_{\tilde{f} \in \mathbb{F}_{\tilde{k}} \cap \partial \mathbb{F}} h_{\tilde{k} \tilde{f}}^{n+1} h_{\tilde{k}}^{n+1} (\overline{w}_{\tilde{k}}^{n\star} - \partial_t^{n+1} B_{\tilde{k}}) \mathcal{N}_{\tilde{k}}^{\tilde{k}_f} |\tilde{f}| \right). $$

and $XX = GN$ for $(GN^\delta)$, the parameters defined in $(GN^\delta.b)$ and

$$ \mu_{\tilde{k}}^{\text{GN}} = \mu_{\tilde{k}}^{\text{NH}} \tilde{f}, \quad \kappa_{\tilde{k}}^{\text{GN}} = \frac{4}{3} \kappa_{\tilde{k}}^{\text{NH}} \tilde{f}, $$

$$ \tilde{\beta}_{\tilde{k}}^{\text{GN}} = \tilde{\beta}_{\tilde{k}}^{\text{NH}} + \frac{1}{4\sqrt{3}|\tilde{k}|} \left( \sum_{f \in \mathbb{F}_{\tilde{k}} \setminus \partial \mathbb{F}} \left( h_{\tilde{k}_f}^{n+1} \right)^2 \sigma_{\tilde{k}_f}^{n\star} \mathcal{N}_{\tilde{k}_f} |f| + \sum_{\tilde{f} \in \mathbb{F}_{\tilde{k}} \cap \partial \mathbb{F}} h_{\tilde{k} \tilde{f}}^{n+1} h_{\tilde{k}}^{n+1} \sigma_{\tilde{k}_f}^{n\star} \mathcal{N}_{\tilde{k}_f}^{\tilde{k}_f} |\tilde{f}| \right). $$

## 3.5 Strategies for improving efficiency

### 3.5.1 Second-order extension

First order schemes are generally not efficient enough for applications. To obtain sufficiently accurate results, it requires a large number of control volumes and therefore a prohibitive computing time. High order strategies can be obtained from the first order scheme, using for example the discontinuous-Galerkin method or the Muscl method for space and the Range-Kutta scheme in time. Note that they already exist high order scheme for the resolution of the dispersive models. For the Green-Naghdi model (GN), we can refer to [30] for Finite Element based method, to [14, 29] for Discontinuous-Galerkin based method and to [34] for residual based method. Here a “second order” extension is proposed by solving the advection step using a Muscl method with flux limiter, and a Runge-Kutta 2 time scheme [40] and a second order hydrostatic reconstruction [3]. The dispersive step is actually already second order in space because of the centered discrete operator (3). In addition, it acts as a projection, see Proposition 5 and Proposition 8. In order to limite the
number of linear system solved by time step, the incremental pressure-correction strategy is used, see [23, 39]. Note no proof of consistency's order is provided in the current work. The “second order” scheme is expected second ordre in $L^2$-norm (time and space) by using classical second order tools. Numerical experiments confirm the expected second order convergence rate, see §4.1.

In the following, we describe more precisely the “second order” scheme. A Runge-Kutta 2 time scheme with Muscl reconstruction for space discretization are used to estimate the solution neglecting the constrains (NH.b) and (GN.b). Note that any second order strategy can be a priori used at this step. The dispersive source terms are added using the estimation of the hydrodynamic pressure at the previous time step. In our case we write with $XX \in \{\text{NH}, \text{GN}\}$

$$U_k^{XX,n+1/2} = U_k^{XX,n} - \frac{\delta_t^n}{2|k|} \sum_{f \in \mathcal{F}_k} \left( \mathcal{F}_f^{XX,n} \cdot \mathcal{N}_k^{k_f} + \mathcal{L}_{f,k}^{XX,n} \right) |f| + \frac{\delta_t^n}{2} \mathcal{D}_k^{XX,n+1/2} (\bar{q}^n, q_B^n)$$

$$U_k^{XX,n*} = U_k^{XX,n} - \frac{\delta_t^n}{|k|} \sum_{f \in \mathcal{F}_k} \left( \mathcal{F}_f^{XX,n+1/2} \cdot \mathcal{N}_k^{k_f} + \mathcal{L}_{f,k}^{XX,n+1/2} \right) |f| + \delta_t^n \mathcal{D}_k^{XX,n+1} (\bar{q}^n, q_B^n)$$

with second order numerical flux $\mathcal{F}_f^{XX,n} = \mathcal{F}^{XX} (U_k^{XX,n}, U_{\bar{k}_f}^{XX,n})$ and source term $\mathcal{L}_{f,k}^{XX,n} = \mathcal{L}^{XX} (U_{\bar{k}_f}^{XX,n}, U_k^{XX,n})$ are computed using the reconstructed states at the face and with a flux limiter (here we use the van Leer limiter), see [3] for more details, and the dispersive operator $\mathcal{D}_k^{XX,n}$ are given in (13) and (23).

Then the solution at time iteration $n+1$ is computed using (NH$^\delta$.b) or (GN$^\delta$.b). The pressure at the next time step is computed considering the dispersion step, i.e.

$$q_{B,k}^{n+1} = q_{B,k}^{n} + \tilde{q}_{B,k}^{n+1} \quad \text{with} \quad \tilde{q}_{B,k}^{n+1} = h_k^{n+1} \frac{\bar{w}_k^{n+1} - \bar{w}_k^{n*}}{\delta_t^n}$$
$$\text{and} \quad \bar{q}_k^{n+1} = \bar{q}_k^n + \tilde{\bar{q}}_k^{n+1} \quad \text{with} \quad \tilde{\bar{q}}_k^{n+1} = \frac{\tilde{q}_{B,k}^{n+1}}{2} + \frac{h_k^{n+1} \sigma_k^{n+1} - \sigma_k^{n*}}{2\sqrt{3} \delta_t^n}$$

It is enough to ensure that for continuous solution the approximation $U^{XX,n*}$ is actually close enough (in order of $(\delta_t^n)^2$) to the set $\mathcal{A}_0^{XX,n+1}$. Therefore the dispersion step (NH$^\delta$.b) or (GN$^\delta$.b), which is actually a projection see Proposition 5 and Proposition 8, introduce an error up to $(\delta_t^n)^2$. See [23] for more details. Eventually, it can be noted that the scheme requires an estimate of the initial pressure, which is unclear. We propose to perform at the initialization one time step with the first order scheme to get an estimation of the hydrodynamic pressure, then to compute again the first time step using the second order scheme.

**Remark 5.** The strategy is particularly efficient since it does not require more resolution of linear system than the first order scheme. Anyway the time of computation for the same mesh is significantly larger since the required time step for the second order shallow water scheme is smaller (generally half as small).
<figure><img src="image_3.png" /><figcaption>Figure 2: §4.1 *Solitary waves* – Water depth approximated by (A.1) and ($GN^\delta$) at time $t = 50$ initialized by the solitary wave (29) of (GN).</figcaption></figure>

### 3.5.2 Adaptive strategy

Because of the implicit resolution, the dispersion step is clearly the more expensive step in term of computational time. Since this step is actually a projection, see Proposition 5 and Proposition 8, it could be relevant to not proceed to the projection step if the solution after the advection step is close enough to the set of admissible solutions. More precisely, let us fix a tolerance $r$ such that instead of looking at a solution exactly satisfying the constrain, i.e. in $A_0^{xx,n}$, we just settle for solution in $A_r^{xx,n}$, see the definition (19) and (28). Then if we have $(\bar{u}^{n*}, \bar{w}^{n*}, \sigma^{n*}) \in A_r^{xx,n+1}$, then we set $(\bar{u}^{n+1}, \bar{w}^{n+1}, \sigma^{n+1}) = (\bar{u}^{n*}, \bar{w}^{n*}, \sigma^{n*})$ and $\bar{q}_k^{n+1} = q_{B,k}^{n+1} = 0$ (or $\tilde{q}_k^{n+1} = \bar{q}_{B,k}^{n+1} = 0$ for the “second order” scheme).

## 4 Numerical results

This section is devoted to numerical illustrations and comparisons between the models in one dimensional framework on a regular grid, i.e. for any $k \in \mathbb{T}$, $|k| = \delta_x$ fixed and defined further. The gravitational acceleration is set to $g = 9.81$ for all following simulations. Unless otherwise specified, the advection step is computed using the HLL scheme [5] with hydrostatic reconstruction [2].

Let us motivate the choice of the test cases. In §4.1, the numerical convergence rates are computed on an analytical solution. In §4.2, the comparison between the entropy-satisfying scheme (GN$^δ$) and the non entropy-satisfying scheme (A.1) is performed. In §4.3, the behavior of the models when the aspect ratio goes to zero is illustrated. In §4.4, a simulation with moving bottom and dry areas is realized to illustrate the robustness of the strategy.

### 4.1 Solitary waves

It is well-known that the dispersive models (NH) and (GN) has a solitary wave
<figure><img src="image_4.png" /><figcaption>Figure 3: §4.1 Solitary waves – L²-error of the numerical approximations (A.1) and (GN^δ), with and without the “second order” strategy, as function of the space step at time t = 50.</figcaption></figure>

solution when $\nabla\phi = 0$. More precisely, the solitary waves of (NH) and (GN) read

$$ (29) \qquad \begin{gathered} h(t,x) = h^0(x) = H_0 + A \operatorname{sech}^2\left(\sqrt{\frac{\gamma^{XX} A}{A+H_0}} \frac{x}{H_0}\right) \\ \bar{u}(t,x) = \bar{u}^0(x) = \frac{H_0}{h^0(x)} \sqrt{g(H_0+A)} \end{gathered} $$

with $H_0 > 0$ is the water depth far from the wave and $A > 0$ is the wave elevation. The parameter $\gamma^{XX}$ depend on the model, i.e. $\gamma^{NH} = 1$ and $\gamma^{GN} = \frac{3}{4}$. The vertical velocity and the hydrodynamic pressure can be deduced using the compatibility conditions (8) and (21).

In the current section, the convergence rate of the schemes (GN^δ) and (A.1) is illustrated in the simulation of the solitary wave (29). More precisely we set the parameters

$$ H_0 = 5 \cdot 10^{-2} \quad \text{and} \quad A = 5 \cdot 10^{-3}. $$

The computational domain is set to [-1, 1]. At the left boundary, the discharged is imposed and set to $h(t, -1) \bar{u}(t, -1) = h^0(-1) \bar{u}^0(-1)$ and at the right the water depth is imposed and set to $h(t, 1) = h^0(1)$.

In Figure 2, the numerical solutions of the schemes, entropic (GN^δ) and compact (A.1), first order and “second order”, are plotted at time $t = 50$ for some space steps. For small enough space step, the solitary wave is well preserved by every schemes. The convergence rate, illustrated in Figure 3, is almost one for all schemes except for the “second order” (GN^δ) scheme for which it is almost two as expected. Note the “second order” strategy does not improve the result of the compact scheme (A.1). This is because the dispersion step (A.1) is not a projection on the set $A_0^{GN,n}$. The same results are observed on the other variables $\bar{u}$, $\bar{w}$ and $\sigma$ and applying the
<figure><img src="image_5.png" /><figcaption>Figure 4: §4.2 Undular jump – Water level approximated (first line) and hydraulic head (second line) for several schemes.</figcaption></figure>

scheme to the non-hydrostatic model (NH$^{δ}$).

## 4.2 Undular jump

A typical test case of (SW) is the transcritical stationary solution in one dimensional framework, see [8]. An analytical solution is based on the energy conservation property, i.e. Hypothesis 2.ii). More precisely, the energy conservation implies that the hydraulic head, defined by the ratio of energy flow to mass flow, i. e.

$$K^{\text{SW}} := \frac{\overline{\mathcal{G}}^h + \overline{\mathcal{G}}^{\overline{u}}}{h\overline{u}} = \phi + gh + \frac{|\overline{u}|^2}{2},$$

is piecewise constant. At the discontinuities, the jump should satisfy the Rankine-Hugoniot relation, see [8] for details. Using the same arguments, the hydraulic head of the dispersive models (NH) and (GN) respectively given by

$$K^{\text{NH}} := \frac{\overline{\mathcal{G}}^h + \overline{\mathcal{G}}^{\overline{u}} + \overline{\mathcal{G}}^{\overline{w}} + \overline{\mathcal{G}}^{q_B}}{h\overline{u}} = \phi + gh + \frac{|\overline{u}|^2}{2} + \frac{|\overline{w}|^2}{2} + \frac{q_B}{2}$$

$$\text{and } K^{\text{GN}} := \frac{\overline{\mathcal{G}}^h + \overline{\mathcal{G}}^{\overline{u}} + \overline{\mathcal{G}}^{\overline{w}} + \overline{\mathcal{G}}^{\sigma} + \overline{\mathcal{G}}^{\overline{q}}}{h\overline{u}} = \phi + gh + \frac{|\overline{u}|^2}{2} + \frac{|\overline{w}|^2}{2} + \frac{|\sigma|^2}{2} + \overline{q}$$

are constant (there is no discontinuity for the dispersive models). Unfortunately, it is not possible to recover the water depth form the hydraulic head as it is for (SW). However, it is possible to compare the hydraulic head of the computed solution to the analytical value. The entropy-stability of the scheme is illustrated by
the monotony of the hydraulic head at steady state.

Let us describe the test case more precisely. We consider the computational domain $[0, 1]$ with an imposed discharged at the left bound $h(t, 0) \bar{u}(t, 0) = 10^{-2}$ and an imposed water depth at the right bound $h(t, 1) = 10^{-1}$. The bottom and the surface pressure are set to

$$B(x) = 7.5 \cdot 10^{-2} e^{-10^2(x-0.3)^2} \quad \text{and} \quad P(x) = 0.$$

Note that the converged solution is out of reach. More precisely, the better the resolution, the more oscillations extend to the right. The converged solution presents periodic oscillations at the right of the obstacle which does not decrease in amplitude. Thus for a fine enough resolution, the right boundary condition (fixed water depth) is no more relevant and more elaborated boundary condition are required (fixed hydraulic head). Unfortunately, it is not clear how to prescribe it at the discrete level.

In Figure 4, the water level and the hydraulic head computed with several schemes $YYY(XX^\delta)$ are plotted for $\delta_x = 10^{-4}$. The schemes $YYY$ refer to the advection step. The CPR (Centered Potential Regularization) scheme [31] is entropy-satisfying, i.e. it satisfies Hypothesis 1.ii), whereas the HLL scheme with hydrostatic reconstruction satisfies only a semi-discrete counterpart of Hypothesis 1.ii). The schemes $(XX^\delta)$ refer to the dispersive step. $(GN^\delta)$ is entropy-satisfying, see Proposition 7.ii), contrary to (A.1). Only the solution CPR$(GN^\delta)$ is fully entropy-satisfying and in practice, it is the most accurate. However, the solution of HLL$(GN^\delta)$ is almost similar. For this both scheme, the hydraulic head is decreasing as a consequence of the mechanical energy dissipation. In the case of (A.1), the hydraulic head does not decrease, the mechanical energy is not well dissipated and the results are significantly less accurate. The impact on the water level is a strong numerical diffusion for both $YYY$ schemes where the hydrodynamic effects occur.

## 4.3 Water drop

The following test case is devoted to the parametric analysis of the aspect ratio. More precisely, (SW) is known as to be a good approximation of (E) when the aspect ratio $\epsilon$ is small enough. We propose to compare the dispersive model (NH) and (GN) to (SW) with respect to $\epsilon$.

Let us describe the case precisely. The computational domain is set to $[0, 1]$ and the boundary conditions are walls (symmetric flow). There is no forcing $\nabla\phi = 0$ and the initial condition is set to

$$h^0(x) = \epsilon (1 + e^{-100x^2}) \quad \text{and} \quad \bar{u}^0(x) = 0.$$

In Figure 5, the water depth approximated by each models with $\delta_x = 10^{-6}$ is plotted for several aspect ratio. The numerical diffusion probably significantly affect the solution for $\epsilon \le 2 \cdot 10^{-4}$. The solution of the dispersive models (NH) and (GN) are very similar, i.e. the solutions present a dispersive shock with oscillations of about
<figure><img src="image_6.png" /><figcaption>Figure 5: §4.3 Water drop – Rescaled water depth $h_\epsilon$ approximated at the rescaled time $\sqrt{g\epsilon t} = 0.6$ by (NH$^\delta$) (left column) and (GN$^\delta$) (right column) for several aspect ratio.</figcaption></figure>

twice the amplitude of the shallow water discontinuity. The frequency is slightly higher for (NH) than for (GN) even if the shape of the convex hull is the same, i.e. linear with almost the same slope. The smaller the $\epsilon$, the higher the frequency of the oscillation. In Figure 6, the CPU time for each simulations is given. The slightly
<figure><img src="image_7.png" /><figcaption>Figure 6: §4.3 Water drop – CPU time of the hierarchy of models to reach the rescaled time $\sqrt{g\epsilon t} = 0.6$ for several aspect ratio.</figcaption></figure>

larger CPU time of (GN$^{δ}$) can be explain by the advection of the standard deviation.

## 4.4 Seawall

To illustrate the robustness of the numerical strategy for moving bottom and with the dry front, the test case of the seawall is presented. The computational domain is set to [0, 10] and the boundary conditions are walls. The surface pressure is neglected $P = 0$ but the bottom is a time and space function given by

$$B(t, x) = \max(0, \min(0.2(x-2), 1)) + 1.5e^{-10*(x-7)^2} + 10e^{-5(x-\min(2t, 2))^2}$$

and the initial condition reads

$$h^0(x) = \begin{cases} \max(0, 2 - B(0, x)) & \text{if } x \le 7 \\ 0 & \text{elsewhere} \end{cases} \quad \text{and} \quad \bar{u}^0(x) = 0.$$

In Figure 7, the water level approximated with $HLL(XX^\delta)$ with $XX \in \{NH, GN\}$ and $\delta_x = 10^{-3}$ is plotted at several times. Even if the two dispersive models lead to significantly different results, especially in the dry front, they are qualitatively similar. First, we find that the dispersive models and (SW) react very differently to the bottom elevation, i.e. $0 \le t < 1$. In the case of (SW), the water level does not become much higher than the initial condition, i.e. 3, and a wave is quickly generated. In the case of dispersive models, the water level becomes higher, i.e. it reaches 4, and the beginning of the wave is not as clear. All models pass the dike, but the shape of the free surface is different. In the case of (SW), the free surface is discontinuous before the dike and continuous after, whereas for dispersive models the free surface is continuous before the dike and discontinuous at the dry front, at least for a while. More precisely, a spike appears at the dry front for a time after the dike and it is clearly not a numerical instability. It was checked that this spike is not very sensitive the velocity in the dry cell. We recall that the velocity in the dry cell is set to zero to impose the
<figure><img src="image_8.png" /><figcaption>Figure 7: §4.4 *Seawall* – Water level approximated by the hierarchy of models at several times.</figcaption></figure>

dissipation of the entropy even at dry front, see Proposition 4 and Proposition 7. It
was checked that for each models, the mechanical energy is decreasing from time
t = 1 (when the bottom stops to move).
## Conclusion

This article focuses on numerical resolution in the multidimensional framework of a hierarchy of vertical-integrated dispersive models. A particular attention is paid to the dissipation of mechanical energy at the discrete level, acting as a mathematical entropy. In terms of efficiency, a numerical improvement of the numerical scheme to obtain a second order convergence rate is proposed using classical tools. In addition, a adaptive strategy is proposed to not solve the implicit part of the scheme when it is not required. To illustrate the accuracy and the robustness of the strategy, several numerical experiments are carried out. In particular, the strategy is capable of treating dry areas.

The first perspective of this work concerns the boundary conditions. If simple boundary conditions can be prescribed as long as the hydrodynamic effects are small at the bound, some other conditions such are transparent condition or incoming wave are also used in practice. The boundary condition that fix the hydraulic head is also required to simulate the hydraulic jump §4.2, which seems a fine test case to compare a model to a complex physics. The second important improvement should be the analysis of the numerical equilibrium of the steady state. The lake at rest is clearly stable assuming that it is stable for $(SW^\delta)$. However, it could be interesting to look for numerical strategy that exactly preserved at the discrete level the solitary wave §4.1. Last but not least, multi-dimensional simulations of the present strategy is also interesting. If the stability property establish in the present work is also valid in horizontal 2D framework, 2D simulations allowed the study of the dispersive models for the shoaling process and the convergence of the waves at the beach. In the vertical direction, a layerwise version of the dispersive models (NH) and (GN) was proposed in [17]. These model can be see as vertical-discretization of the 3D Euler model with free surface (E). The adaptation of the numerical strategy proposed in the current work to the layerwise version of the dispersive models should lead to a entropy-satisfying scheme of (E).

## Acknowledgements

I would like to thank the Alpines Inria project team members with in particular S. Cayrols, O. Tissot, H. Al Daas for their help with the linear algebra library (Intel® Math Kernel Library) used for this work, and E. Audusse for his kind advises. This work is a part of the research project CNRS-PICS-07480.

## A Dispersion step with a compact operator

$(NH^\delta.b)$ or $(GN^\delta.b)$ are a 5-point stencil operator in one dimensional framework. One can remark that the scheme is a discretization of a reaction-advection-diffusion equation. A strategy that can be use to reduce the stencil of the scheme is to establish the reaction-advection-diffusion equation mimicking the step realized to obtain the scheme $(NH^\delta.b)$ or $(GN^\delta.b)$ at the continuous level, then use a more compact
discretization. The resulting scheme reads

$$ (A.1) \qquad \alpha_k^{xx} \bar{u}_k^{n+1} + \nabla_k^\delta (\mu^{xx} \cdot \bar{u}^{n+1}) - \mu_k^{xx} \nabla_k^\delta \bar{u}^{n+1} \\ - \frac{1}{|k|} \sum_{f \in \mathcal{F}_k} \left( \frac{\kappa_k^{xx}}{\delta_k} + \frac{\kappa_{\bar{k}_f}^{xx}}{\delta_{\bar{k}_f}} \right) \frac{\bar{u}_{\bar{k}_f}^{n+1} - \bar{u}_k^{n+1}}{2} \cdot \mathcal{N}_k^{\bar{k}_f} \mathcal{N}_{\bar{k}_f}^{\bar{k}_f} |f| = \beta_k^{xx} $$

with the parameters $\alpha_k^{xx}$, $\mu_k^{xx}$, $\kappa_k^{xx}$ and $\beta_k^{xx}$ ($\mathcal{X} \in \{\text{NH}, \text{GN}\}$) defined in (NH$^\delta$.b) or (GN$^\delta$.b) and $\delta_k$ is the compactness defined in (6). The yielding scheme (A.1) is more compact, i.e. 3-point stencil in one dimension. However, the compact scheme (A.1) does not ensure the mechanical energy dissipation and in less accurate than (NH$^\delta$.b) or (GN$^\delta$.b) in practice, see §4.2.

## References

[1] N. AISSIOUENE, M.-O. BRISTEAU, E. GODLEWSKI, AND J. SAINTE-MARIE, A combined finite volume - finite element scheme for a dispersive shallow water system, Networks and Heterogeneous Media (NHM), (2016).

[2] E. AUDUSSE, F. BOUCHUT, M.-O. BRISTEAU, R. KLEIN, AND B. PERTHAME, A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows, SIAM J. Sci. Comput., 25 (2004), pp. 2050–2065.

[3] E. AUDUSSE AND M.-O. BRISTEAU, A well-balanced positivity preserving "second-order" scheme for shallow water flows on unstructured meshes, Journal of Computational Physics, 206 (2005), pp. 311–333.

[4] P. BONNETON, F. CHAZEL, D. LANNES, F. MARCHE, AND M. TISSIER, A splitting approach for the fully nonlinear and weakly dispersive Green–Naghdi model, Journal of Computational Physics, 230 (2011), pp. 1479–1498.

[5] F. BOUCHUT, Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources, Springer Science & Business Media, 2004.

[6] M.-O. BRISTEAU, A. MANGENEY, J. SAINTE-MARIE, AND N. SEGUIN, An energy-consistent depth-averaged Euler system: Derivation and properties, Discrete and Continuous Dynamical Systems - Series B, 20 (2015), pp. 961–988.

[7] C. CANCES, H. MATHIS, AND N. SEGUIN, Error estimate for time-explicit finite volume approximation of strong solutions to systems of conservation laws, SIAM Journal on Numerical Analysis, 54 (2016), pp. 1263–1287.

[8] H. CHANSON, Open-channel hydraulics, Butterworth-Heinemann, 1959.

[9] F. CHAZEL, D. LANNES, AND F. MARCHE, Numerical Simulation of Strongly Nonlinear and Dispersive Waves Using a Green–Naghdi Model, Journal of Scientific Computing, 48 (2010), pp. 105–116.
[10] R. CIENFUEGOS, E. BARTHÉLEMY, AND P. BONNETON, *A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part I: model development and analysis*, International Journal for Numerical Methods in Fluids, 51 (2006), pp. 1217-1253.

[11] F. COQUEL, J.-M. HÉRARD, K. SALEH, AND N. SEGUIN, *A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model*, ESAIM: Mathematical Modelling and Numerical Analysis, 48 (2013), pp. 165-206.

[12] C. COURTÈS, F. LAGOUTIÈRE, AND F. ROUSSET, *Error estimates of finite difference schemes for the Korteweg-de Vries equation*, (2018).

[13] A.-J.-C. B. DE SAINT-VENANT, *Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leurs lits*, C.R. Acad. Sci. Paris, 73 (1871), pp. 147-154.

[14] A. DURAN AND F. MARCHE, *Discontinuous-galerkin discretization of a new class of green-naghdi equations*, Communications in Computational Physics, 17 (2015), pp. 721-760.

[15] R. DUTTA, U. KOLEY, AND N. RISEBRO, *Convergence of a higher order scheme for the korteweg-de vries equation*, SIAM Journal on Numerical Analysis, 53 (2015), pp. 1963-1983.

[16] N. FAVRIE AND S. GAVRILYUK, *A rapid numerical method for solving serre-green-naghdi equations describing long free surface gravity waves*, Nonlinearity, 30 (2017), p. 2718.

[17] E. D. FERNÁNDEZ-NIETO, M. PARISOT, Y. PENEL, AND J. SAINTE-MARIE, *A hierarchy of dispersive layer-averaged approximations of Euler equations for free surface flows*, Communications in Mathematical Sciences, 16 (2018), pp. 1169-1202.

[18] T. GALLOUËT, J.-M. HÉRARD, AND N. SEGUIN, *Numerical modeling of two-phase flows using the two-fluid two-pressure approach*, Mathematical Models and Methods in Applied Sciences, 14 (2004), pp. 663-700.

[19] J.-F. GERBEAU AND B. PERTHAME, *Derivation of viscous saint-venant system for laminar shallow water; numerical validation*, Discrete and Continuous Dynamical Systems - Series B, 1 (2001), pp. 89-102.

[20] E. GODLEWSKI AND P.-A. RAVIART, *Numerical approximation of hyperbolic systems of conservation laws*, vol. 118 of Applied Mathematical Sciences, Springer-Verlag, New York, 1996.

[21] S. K. GODUNOV, *A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics*, Mat. Sb. (N.S.), 47(89) (1959), pp. 271-306.
[22] A. E. GREEN AND P. M. NAGHDI, *A derivation of equations for wave propagation in water of variable depth*, Journal of Fluid Mechanics, 78 (1976), pp. 237-246.

[23] J. GUERMOND, P. MINEV, AND J. SHEN, *An overview of projection methods for incompressible flows*, Computer Methods in Applied Mechanics and Engineering, 195 (2006), pp. 6011–6045.

[24] H. HOLDEN, C. LUBICH, AND N. H. RISEBRO, *Operator splitting for partial differential equations with burgers nonlinearity*, Mathematics of Computation, 82 (2012), pp. 173–185.

[25] D. LANNES AND P. BONNETON, *Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation*, Physics of Fluids, 21 (2009).

[26] D. LANNES AND F. MARCHE, *A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations*, Journal of Computational Physics, 282 (2015), pp. 238-268.

[27] O. LE MÉTAYER, S. GAVRILYUK, AND S. HANK, *A numerical scheme for the Green-Naghdi model*, Journal of Computational Physics, 229 (2010), pp. 2034–2045.

[28] R. J. LEVEQUE, *Finite volume methods for hyperbolic problems*, vol. 31, Cambridge university press, 2002.

[29] M. LI, P. GUYENNE, F. LI, AND L. XU, *High order well-balanced cdg-fe methods for shallow water waves by a green-naghdi model*, Journal of Computational Physics, 257 (2014), pp. 169–192.

[30] D. MITSOTAKIS, B. ILAN, AND D. DUTYKH, *On the galerkin/finite-element method for the serre equations*, Journal of Scientific Computing, 61 (2014), pp. 166–195.

[31] M. PARISOT AND J.-P. VILA, *Centered-potential regularization for the advection upstream splitting method*, SIAM Journal on Numerical Analysis, 54 (2016), pp. 3083–3104.

[32] D. H. PEREGRINE, *Long waves on a beach*, Journal of Fluid Mechanics, 27 (1967), pp. 815-827.

[33] B. PERTHAME AND C. SIMEONI, *A kinetic scheme for the saint-venant system with a source term*, Calcolo, 38 (2001), pp. 201-231.

[34] M. RICCHIUTO AND A. FILIPPINI, *Upwind residual discretization of enhanced boussinesq equations for wave propagation over complex bathymetries*, Journal of Computational Physics, 271 (2014), pp. 306–341. Frontiers in Computational Physics.

[35] F. SERRE, *Contribution à l'étude des écoulements permanents et variables dans les canaux*, La Houille Blanche, (1953), pp. 830-872.
[36] J. A. SETHIAN AND P. SMEREKA, *Level set methods for fluid interfaces*, Annual Review of Fluid Mechanics, 35 (2003), pp. 341–372.

[37] I. SULICIU, *On modelling phase transitions by means of rate-type constitutive equations. shock wave structure*, International Journal of Engineering Science, 28 (1990), pp. 829–841.

[38] E. F. TORO, *Riemann solvers and numerical methods for fluid dynamics: a practical introduction*, Springer Science & Business Media, 2013.

[39] J. VAN KAN, *A second-order accurate pressure-correction scheme for viscous incompressible flow*, SIAM Journal on Scientific and Statistical Computing, 7 (1986), pp. 870–891.

[40] B. VAN LEER, *Towards the ultimate conservative difference scheme. v. a second-order sequel to godunov's method*, Journal of Computational Physics, 32 (1979), pp. 101–136.

[41] Y. YAMAZAKI, Z. KOWALIK, AND K. F. CHEUNG, *Depth-integrated, non-hydrostatic model for wave breaking and run-up*, International Journal for Numerical Methods in Fluids, 61 (2009), pp. 473–497.