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

HAL
open science

# Numerical simulations of a dispersive model approximating free-surface Euler equations

Cipriano Escalante-Sánchez, Enrique Domingo Fernandez-Nieto, Tomas

Morales de Luna, Yohan Penel, Jacques Sainte-Marie

► To cite this version:

Cipriano Escalante-Sánchez, Enrique Domingo Fernandez-Nieto, Tomas Morales de Luna, Yohan Penel, Jacques Sainte-Marie. Numerical simulations of a dispersive model approximating free-surface Euler equations. 2020. hal-03079571

HAL Id: hal-03079571

https://hal.science/hal-03079571

Preprint submitted on 17 Dec 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Numerical simulations of a dispersive model
approximating free-surface Euler equations

Cipriano Escalante-Sánchez *, Enrique D. Fernández-Nieto, † T. Morales de Luna, ‡
Yohan Penel, § Jacques Sainte-Marie ¶

December 17, 2020

Abstract

In some configurations, dispersion effects must be taken into account to improve the simulation of complex fluid flows. A family of free-surface dispersive models has been derived in [24]. The hierarchy of models is based on a Galerkin approach and parametrised by the number of discrete layers along the vertical axis. In this paper we propose some numerical schemes designed for these models in a 1D open channel. The cornerstone of this family of models is the *Serre – Green-Naghdi* model which has been extensively studied in the literature from both theoretical and numerical points of view. More precisely, the goal is to propose a numerical method for $LDNH_2$ model that can be defined in terms of a projection method for the one-layer case, despite the number of layers. To do so, the monolayer case is addressed by means of a projection-correction method applied to a non-standard differential operator. A special attention is paid to boundary conditions. This case is extended to several layers thanks to an original relabelling of the unknowns. In the numerical tests we show the converge of the method and its accuracy in comparison with $LDNH_0$ model.

<footnote>*Dpto. de Matemáticas. Universidad de Córdoba. Campus de Rabanales. 14071 Córdoba, Spain. (cescalante@uco.es)</footnote>

<footnote>†Dpto. Matemática Aplicada I. ETS Arquitectura - Universidad de Sevilla. Avda. Reina Mercedes N. 2. 41012-Sevilla, Spain.
(edofer@us.es)</footnote>

<footnote>‡Dpto. de Matemáticas. Universidad de Córdoba. Campus de Rabanales. 14071 Córdoba, Spain. (tomas.morales@uco.es)</footnote>

<footnote>§CEREMA, Sorbonne Université, Université Paris-Diderot Sorbonne Paris Cité, CNRS, INRIA, Laboratoire Jacques-Louis Lions, Paris, France (yohan.penel@inria.fr)</footnote>

<footnote>¶Inria Paris, 2 rue Simone If, CS 42112, 75589 Paris Cedex 12, France & Sorbonne Université, Université de Paris, CNRS, Laboratoire Jacques-Louis Lions, LJLL, F-75005 Paris, France (Jacques.Sainte-Marie@inria.fr)</footnote>
# Contents

1 Introduction 2

2 Analysis and numerical method for the LDNH$_2$(1) model 5

2.1 Properties of the model at the continuous level 6

2.2 Splitting strategy at the semi-discrete level 8

2.3 Numerical method at the discrete level 10

2.3.1 Hyperbolic step 10

2.3.2 Non-hydrostatic step: numerical scheme for the velocity-pressure problem 10

2.3.3 Numerical scheme for the pressure problem 13

3 Numerical scheme for the LDNH$_2$(L) model 14

3.1 Notations and model 14

3.2 Differential operators 15

3.3 Structure of the fully discretised equations 16

3.4 Iterative scheme 18

4 Numerical simulations 19

4.1 Convergence test 19

4.2 Solitary wave propagation over reefs 20

4.3 Wave propagation over a submerged bar 20

4.4 Shoaling of a solitary wave on a plane beach 22

5 Conclusion 22

6 Acknowledgements 27

## 1 Introduction

Water waves and more generally water flows are of great interest in several scientific fields with applications to society issues such as protection of populations (tsunamis, floods, ...) or energy production (water-turbines, ...). Depending on the accuracy that is required in the applications, more or less complex models are used, from fully resolved to averaged equations.

In order to simulate the behaviour of free-surface fluid flows, let us consider the 2D Euler equations for an incompressible free-surface flow under gravity:

$$ \left\{ \begin{array}{l} \nabla \cdot \hat{\mathbf{u}} = 0, \\ \partial_t \hat{\mathbf{u}} + (\hat{\mathbf{u}} \cdot \nabla) \hat{\mathbf{u}} + \nabla \hat{p} = (0, -g), \end{array} \right. \qquad (1a) $$

$$ \partial_t \hat{u} + (\hat{u} \cdot \nabla) \hat{u} + \nabla \hat{p} = (0, -g), \qquad (1b) $$

set in the moving domain (see Figure 1)

$$ \Omega(t) = \{(x, z) \in \mathbb{R}^2 \mid x \in I, z_b(x) \le z \le \hat{\eta}(t, x)\}. $$

Here, $I = (x_l, x_r)$ is a bounded interval of $\mathbb{R}$ and $z_b$ is the given topography (independent from time). The unknowns are the velocity field $\hat{\mathbf{u}} = (\hat{u}, \hat{w})$, the pressure field $\hat{p}$ and the water elevation $\hat{\eta}$. The water height is deduced from $\hat{h}(t, x) = \hat{\eta}(t, x) - z_b(x)$. Viscosity effects are not taken into account in this work but we refer to [9] for instance to deal with such terms.
<figure><img src="image_2.png" /><figcaption>Figure 1: Fluid domain for the Euler equations (1)</figcaption></figure>

The model is supplemented with kinematic boundary conditions:

$$ \begin{cases} \partial_t \hat{\eta} + \hat{u}_s \partial_x \hat{\eta} - \hat{w}_s = 0, \\ \hat{u}_b \partial_x z_b - \hat{w}_b = 0, \end{cases} \qquad (2) $$

as well as

$$ \hat{p}(t, x, \hat{\eta}(t, x)) = p^{\text{atm}}(t, x), \qquad (3) $$

for some given atmospheric pressure $p^{\text{atm}}(t, x) > 0$. Classically, the pressure field is decomposed into three parts:

$$ \hat{p}(t, x, z) = g(\hat{\eta}(t, x) - z) + p^{\text{atm}}(t, x) + \hat{q}(t, x, z) $$

where $\hat{q}$ is the hydrodynamic pressure field – or commonly referred to as the non-hydrostatic component. The hydrostatic part is $g(\hat{\eta}(t, x) - z)$. Hydrostatic models such as the nonlinear shallow water equations [14] or the hydrostatic Navier-Stokes equations [3] are based on the assumption $\hat{q} \equiv 0$.

Given this decomposition, BC (3) is equivalent to

$$ \hat{q}(t, x, \hat{\eta}(t, x)) = 0. \qquad (4) $$

In addition to its complex mathematical structure (see e.g. [7, 27, 37]), Model (1) coupled to (2-4) is a real challenge owing to the fact the domain is moving: the water elevation is an unknown in itself. That is why, in spite of the increase of computer performance, reduced-complexity models have been introduced, analysed and discretised. There exists an extensive literature about models approximating the Euler equations under simplifying hypotheses. Depending on the physical phenomena at stake and thus the fluid regime under study, some models turn out to be more accurate than others.

For example, the nonlinear shallow water equations (NLSW) [14, 25, 26, 41] provide relevant results for large wavelengths but seem restrictive in other regimes, in particular due to the absence of dispersive effects. To go further, models were derived in a given regime of magnitude for the nonlinearity ($\varepsilon$: wave amplitude/water depth ratio) and for the frequency dispersion ($\mu$: water depth/wavelength ratio). The competition between the two phenomena is responsible for the shape of water waves. Dispersion is necessary for instance for stratified flows or close to coastal areas. Different regimes lead to weakly/fully nonlinear – weakly/fully dispersive models. For example, the NLSW equations correspond to $\mu = 0$ for any $\varepsilon = \mathcal{O}(1)$.

We may refer to pioneering works e.g. Boussinesq [8] and Peregrine [44], models obtained to improve linear dispersion properties e.g. Madsen & Sorensen [40] and Nwogu [42], models suitable for deep waters [34], ... For reviews of dispersive models, see [31, 33]. A focus will be made below on two systems: the Depth-Averaged Euler (DAE) system (5) described in [1] and the Serre – Green-Naghdi equations (6) (see [28, 47]).
It is worth noticing that classically dispersive models, usually written under a Boussinesq form, introduce high order derivatives for two unknowns, namely water height and velocity). Contrary to them, non-hydrostatic pressure models introduce such effects by means of the non-hydrostatic component of the pressure. The latter is the strategy used here. Classical dispersive systems, and in particular (5) and (6), may be written under the form of non-hydrostatic system. We refer the reader to [20] for the link between the two approaches. Non-hydostatic formulation has several advantages from the numerical point of view, especially regarding boundary conditions. Notice that they do not rely on any irrotational assumption. Moreover, both (5) and (6) are proven in [24] to belong to the same hierarchy of models depending on some orders of approximation for the set of unknowns<fnref n="1" />. Indeed, the key point is the dependence of the solution on the vertical coordinate z.

More precisely, Equations (5) are referred to as $LDNH_0(1)$ as an approximation of the Euler equations where $(u, w, q)$ is a $(\mathbb{P}_0, \mathbb{P}_0, \mathbb{P}_1)$ Galerkin approximation of $(\hat{u}, \hat{w}, \hat{q})$. Likewise, Equations (6) are referred to as $LDNH_2(1)$ and $(u, w, q)$ is a $(\mathbb{P}_0, \mathbb{P}_1, \mathbb{P}_2)$ Galerkin approximation of $(\hat{u}, \hat{w}, \hat{q})$. In particular, (6) has a more accurate linear dispersion relation than (5) for large wavelengths [24, Fig. 2].

To go further into details, when the frequency dispersion is small, the Euler equations can be approximated by the following non-hydrostatic model [1] named DAE

$$ egin{cases} 
\partial_t h + \partial_x(hu) = 0, \\
\partial_t(hu) + \partial_x(hu^2 + hq) = -\frac{\gamma^2}{2}q\partial_x z_b - gh\partial_x\eta - h\partial_x p^{atm}, \\
\partial_t(hw) + \partial_x(huw) = \gamma q, \\
\gamma w + h\partial_x u - \frac{\gamma^2}{2}u\partial_x z_b = 0,
\end{cases} 
\quad (5) $$

where $(u, w)$ is an approximation of $(\hat{u}, \hat{w})$ along the water column, $q$ is the non-hydrostatic component of the pressure field and $\gamma > 0$ a parameter. The case $\gamma = 2$ was derived in [10] and simulated in [2] and corresponds to $LDNH_0(1)$ in [24]. For a flat topography $\partial_x z_b \equiv 0$, $\gamma = \sqrt{3}$ corresponds to the Serre-Green-Naghdi equations.

Another shallow water model is investigated in the present work. More precisely, we consider the dispersive model:

$$ \partial_t h + \partial_x(hu) = 0, \quad (6a) $$

$$ \partial_t(hu) + \partial_x(hu^2) + \partial_x(hq) + q_b \partial_x z_b = -gh\partial_x(z_b + h) - h\partial_x p^{atm}, \quad (6b) $$

$$ \partial_t(hw) + \partial_x(huw) - q_b = 0, \quad (6c) $$

$$ \partial_t(h\sigma) + \partial_x(h\sigma u) - 2\sqrt{3}\left(q - \frac{q_b}{2}\right) = 0, \quad (6d) $$

$$ w - u\partial_x z_b - \sqrt{3}\sigma = 0, \quad (6e) $$

$$ 2\sqrt{3}\sigma + h\partial_x u = 0. \quad (6f) $$

System (6) is made of first-order PDEs which account for conservation of mass (6a) and of momentum (6b-6d). (6e-6f) are diagnostic equations deduced from the incompressibility constraint (1a). There are several equivalent ways to write down the diagnostic equations (6e) and (6f) but this definition is the only choice that provides the duality relation (9) between the modified "pressure gradient" and the modified "velocity divergence" (8b).

This duality property is crucial from both theoretical and numerical points of view as highlighted in the following. It also allows to specify the boundary conditions:

1. From the weak formulation associated to the underlying elliptic equation, we deduce the boundary conditions for the pressure fields (see Remark 3);

2. We then infer from the numerical point of view how boundary conditions should be naturally imposed by symmetry in the resulting linear system (see § 2.3.2).

<footnote n="1">¹The models of the hierarchy are denoted $LDNH_k(L)$ where $LDNH$ stands for Layerwise Discretisation Non-Hydrostatic, $k$ is the order of approximation and $L$ the number of layers.</footnote>
For the sake of simplicity and without restrictions<fnref n="2" />, we will consider in this paper a subcritical flow (which corresponds to the hyperbolic regime<fnref n="3" /> in the Euler system (1)). Hence we prescribe the discharge on the left ($x = x_l$) – see (7c) and the water height on the right ($x = x_r$) – see (7d). Then, Model (6) is supplemented with the following boundary conditions:

* For the non-hydrostatic pressure, we impose:

$$ \text{At } x = x_{\ell}, \qquad \partial_x(hq) = \chi_{\ell}; \qquad (7a) $$

$$ \text{At } x = x_r, \qquad q = 0, \qquad (7b) $$

for some flux $\chi_{\ell} \in \mathbb{R}$. The choice of these boundary conditions relies on the weak formulation associated to the underlying elliptic problem (see Remark 3) and on the boundary condition for the topography (7e).

* For the velocity field and the water height, it depends on the underlying hyperbolic regime. As mentioned above, we consider a subcritical flow so that, following the investigations of § 2.3.2, we impose

$$ \text{At } x = x_{\ell}, \qquad hu = Q_u, \qquad hw = Q_w, \qquad h\sigma = \frac{Q_w}{\sqrt{3}}; \qquad (7c) $$

$$ \text{At } x = x_r, \qquad h = h_r; \qquad (7d) $$

for some constants $Q_u > 0, Q_w \ge 0$ and $h_r > 0$.

We also assume that the topography satisfies the following boundary condition:

$$ \partial_x z_b(x_{\ell}) = 0. \qquad (7e) $$

It is shown in § 2.1 that System (6) can be rewritten under three different other formulations (Props. 1, 2 and 3). In particular, Proposition 1 shows that (6) is equivalent to the *Serre – Green-Naghdi* equations under Boussinesq form, as presented in [33]. Then, from now on, System (6) coupled to boundary conditions (7) will be referred to indifferently as SGN or $LDNH_2(1)$.

Although Model $LDNH_2(1)$ has a larger number of unknowns than $LDNH_0(1)$ (which may provide a richer modelling), the first goal of this paper is to design a numerical method for (6) with a reduced additional computational complexity compared to (5). See [43] for another strategy. The second goal of this paper is to extend the aforementioned numerical method designed for $LDNH_2(1)$ to its multilayer counterpart $LDNH_2(L)$ (for some number of layers $L > 1$) as derived in [24]. Indeed, $LDNH_2(1)$ relies on the approximation that the horizontal component of the velocity field does not depend on $z$, i.e. it is constant along the water column. In some cases (for instance when the flow is not shallow), it is necessary to add more degrees of freedom. The Galerkin method used for the semi-discretisation in $z$ leads to $LDNH_2(L)$ where $L$ is the number of vertical cells.

The paper is organised as follows: in Section 2, we first show some properties of the $LDNH_2(1)$ model: equivalent formulations (including the relationship with the *Serre – Green-Naghdi* equations) and associated energy. Secondly, a numerical strategy to solve the monolayer model $LDNH_2(1)$ is presented: it consists in an iterative algorithm taking into account the gradient-divergence duality. In Section 3, the layer-averaged extension $LDNH_2(L)$ is reformulated with similar differential operators. The previous numerical method is then extended to the multilayer case, still based on the gradient-divergence duality. In Section 4, some classic numerical tests are presented to assess these strategies for $LDNH_2(1)$ and $LDNH_2(L)$, in particular in terms of accuracy with respect to $LDNH_0(L)$.

## 2 Analysis and numerical method for the $LDNH_2(1)$ model

In this section we first study some of the properties of Model (6), such as some reformulations and its associated energy. Secondly, the design of a numerical discretisation of the model is presented.

<footnote n="2">²The location of the boundary condition upon the water height does not interfere with the projection step.</footnote>

<footnote n="3">³The hyperbolic part is not the core of this paper which is rather the elliptic part.</footnote>
## 2.1 Properties of the model at the continuous level

Let us set

$$ \mathbf{X} = \begin{pmatrix} u \\ w \\ \sigma \end{pmatrix}, \quad \mathbf{Q} = \begin{pmatrix} q \\ q_b \end{pmatrix}, \quad \text{and} \quad \mathbf{S}_h = \begin{pmatrix} -gh\partial_x(z_b+h) - h\partial_x p^{\text{atm}}, & 0 \\ 0 & 0 \end{pmatrix}. \qquad (8a) $$

Let us also define the operators

$$ \nabla_{\text{sgn}} \mathbf{Q} = \begin{pmatrix} \partial_x(hq) + q_b \partial_x z_b \\ -q_b \\ -2\sqrt{3}\left(q - \frac{q_b}{2}\right) \end{pmatrix} \quad \text{and} \quad \nabla_{\text{sgn}} \mathbf{X} = \begin{pmatrix} 2\sqrt{3}\sigma + h\partial_x u \\ w - u\partial_x z_b - \sqrt{3}\sigma \end{pmatrix} \qquad (8b) $$

and notice that the following duality relation holds

$$ \mathbf{X} \cdot \nabla_{\text{sgn}} \mathbf{Q} = \partial_x(hqu) - \mathbf{Q} \cdot (\nabla_{\text{sgn}} \mathbf{X}). \qquad (9) $$

Given these notations, System (6) reads in a compact form

$$ \left\{ \begin{aligned} & \partial_t h + \partial_x(hu) = 0, \\ & \partial_t(h\mathbf{X}) + \partial_x(hu\mathbf{X}) + \nabla_{\text{sgn}} \mathbf{Q} = \mathbf{S}_h, \\ & \nabla_{\text{sgn}} \mathbf{X} = 0. \end{aligned} \right. \qquad (10a) \\ (10b) \\ (10c) $$

System (10) has the same mathematical structure as the incompressible Euler equations with variable density except that differential operators are non-standard. Hence, this remark will help to design a similar numerical strategy.

It is shown below that this system can be rewritten under three different other forms (Props. 1, 2 and 3).

First, let us remark that System (10) is nothing but the Serre-Green-Naghdi equations as presented in [33]:

**Proposition 1.** System (10) can be rewritten under the Boussinesq formulation

$$ \left\{ \begin{aligned} & \partial_t h + \partial_x(hu) = 0, \\ & (\mathcal{I}_d + \mathcal{T}[h, z_b])(\partial_t u + u\partial_x u) + g\partial_x(h + z_b) + \mathcal{Q}[h, z_b]u = -\partial_x p^{\text{atm}}, \end{aligned} \right. \qquad (11) $$

where

$$ \mathcal{T}[h, z_b]v = \mathcal{R}_1[h, z_b](\partial_x v) + \mathcal{R}_2[h, z_b](v \partial_x z_b), $$

$$ \mathcal{Q}[h, z_b]v = -2\mathcal{R}_1[h, z_b]((\partial_x v)^2) + \mathcal{R}_2[h, z_b](v^2 \partial_{xx}^2 z_b), $$

$$ \mathcal{R}_1[h, z_b]w = -\frac{1}{3h}\partial_x(h^3 w) - \frac{h}{2}w \partial_x z_b, $$

$$ \mathcal{R}_2[h, z_b]w = \frac{1}{2h}\partial_x(h^2 w) + w \partial_x z_b. $$

*Proof.* Indeed, we deduce from (6) that

$$ \sigma = -\frac{h \partial_x u}{2\sqrt{3}}, \quad w = -\frac{h \partial_x u}{2} + u \partial_x z_b, $$

$$ q_b = -\frac{h^2}{2} \left[ \partial_x (\partial_t u + u \partial_x u) - 2(\partial_x u)^2 \right] + h \left[ \partial_x z_b (\partial_t u + u \partial_x u) + u^2 \partial_{xx}^2 z_b \right], $$

$$ q = -\frac{h^2}{3} \left[ \partial_x (\partial_t u + u \partial_x u) - 2(\partial_x u)^2 \right] + \frac{h}{2} \left[ \partial_x z_b (\partial_t u + u \partial_x u) + u^2 \partial_{xx}^2 z_b \right]. $$

Inserting the two last equalities into (6b) leads to the expected result. □
The main difference between the non-hydrostatic formulation (6) and the Boussinesq formulation given in Proposition 1 lies in the order of spatial derivatives (1 vs. 3). The consequence is a larger number of unknowns in (6) (6 vs. 2) but with at most first order derivatives. One key advantage of (6) is that the treatment of boundary conditions is easier as it will be shown later. Moreover, a smaller stencil is needed when one has to approach numerically first order compared to higher order derivatives.

A numerical algorithm to simulate (11) is designed in [6] by introducing a second-order parametrised perturbation and based on the inversion of $\mathcal{I}_d + \mathcal{T}$. The parameter is set so that the linear dispersion relation is optimised with respect to the Airy relation. The numerical technique consists of a splitting method between the hyperbolic Saint-Venant equations (solved with a finite-volume scheme) and the dispersive part (solved with a finite-difference scheme). An extension to dimension 2 is provided in [35] by replacing the differential operator $\mathcal{I}_d + \mathcal{T}$ by a time independent “diagonal” approximation. In [16], a Discontinuous Galerkin Finite-Element method is applied.

We also mention the following formulation with high order derivatives of $h$:

**Proposition 2.** System (10) also reads

$$ \left\{ \begin{aligned} \partial_t h + \partial_x(hu) &= 0, \\ \partial_t(hu) + \partial_x \left( hu^2 + g \frac{h^2}{2} + \frac{h^2 \ddot{h}}{3} + \frac{h^2 \dot{u}}{2} \partial_x z_b + \frac{h^2 u^2}{2} \partial_{xx}^2 z_b \right) \\ &\quad + h \left( g + \frac{\ddot{h}}{2} + \dot{u} \partial_x z_b + u^2 \partial_{xx}^2 z_b \right) \partial_x z_b &= -h \partial_x p^{\text{atm}}. \end{aligned} \right. $$

Here we used the standard notation $\dot{\xi} := \partial_t \xi + u \partial_x \xi$.

*Proof.* We deduce from (6)

$$ \begin{aligned} \dot{h} &= -h \partial_x u = 2\sqrt{3}\sigma = 2(w - u \partial_x z_b), \\ q_b &= h \dot{w} = \frac{h \ddot{h}}{2} + h \left[ \partial_t(u \partial_x z_b) + u \partial_x(u \partial_x z_b) \right], \\ q &= \frac{q_b}{2} + \frac{h \ddot{h}}{12} = \frac{h \ddot{h}}{3} + \frac{h}{2} \left[ \partial_t(u \partial_x z_b) + u \partial_x(u \partial_x z_b) \right]. \end{aligned} $$

Inserting the two latter equalities into (6b), we obtain the expected result. □

**Remark 1.** For a flat topography, we recover the model studied in [22,36]. In [36], the model is reformulated by means of the change of variable $(h, hu) \mapsto (h, u + \frac{\partial_x(h^2 \dot{h})}{3h})$ and solved using a finite-difference method. The numerical technique used in [22] is based on an augmented-Lagrangian approach.

In the framework of projection methods, we exhibit another formulation. The two constraints (10c) can be replaced by applying the divergence (8b) to the momentum equation (10b) (initially divided by $h$):

**Proposition 3.** System (10) can be rewritten under the reduced formulation

$$ \partial_t h + \partial_x(hu) = 0, \qquad (12a) $$

$$ \partial_t(h\mathbf{X}) + \partial_x(hu\mathbf{X}) + \nabla_{\text{sgn}} \mathbf{Q} = \mathbf{S}_h, \qquad (12b) $$

$$ - \nabla_{\text{sgn}} \cdot \left( \frac{\nabla_{\text{sgn}} \mathbf{Q}}{h} \right) = - \nabla_{\text{sgn}} \cdot \left( \frac{\mathbf{S}_h}{h} \right) + \nabla_{\text{sgn}} \cdot (u \partial_x \mathbf{X}). \qquad (12c) $$

Eq. (12c) reads

$$ 12 \frac{q}{h} - h \partial_x \left( \frac{\partial_x(hq)}{h} \right) - 6 \frac{q_b}{h} - h \partial_x \left( \frac{q_b}{h} \partial_x z_b \right) = 2h(\partial_x u)^2 + h \partial_x (g \partial_x(z_b + h) + \partial_x p^{\text{atm}}), \qquad (13a) $$

$$ (4 + (\partial_x z_b)^2) \frac{q_b}{h} - 6 \frac{q}{h} + \partial_x z_b \frac{\partial_x(hq)}{h} = u^2 \partial_{xx}^2 z_b - (g \partial_x(z_b + h) + \partial_x p^{\text{atm}}) \partial_x z_b. \qquad (13b) $$
Well-posedness of this system is studied in [4].

**Remark 2.** Notice that only h and u are involved (not w nor σ) in the right hand side of (13). Moreover, we also notice that $q_b$ can be expressed directly from (13b)

$$ \frac{q_b}{h} = \frac{1}{4 + (\partial_x z_b)^2} \left[ 6 \frac{q}{h} - \partial_x z_b \frac{\partial_x(hq)}{h} + u^2 \partial_{xx}^2 z_b - (g \partial_x(z_b + h) + \partial_x p^{atm}) \partial_x z_b \right]. \quad (14) $$

Inserting the latter expression into (13a) provides a unique equation for q:

$$ 12 \frac{1 + (\partial_x z_b)^2}{4 + (\partial_x z_b)^2} \frac{q}{h} - h \partial_x \left( \frac{4}{4 + (\partial_x z_b)^2} \frac{\partial_x(hq)}{h} \right) + \frac{6 \partial_x z_b}{4 + (\partial_x z_b)^2} \frac{\partial_x(hq)}{h} - h \partial_x \left( \frac{6 \partial_x z_b}{4 + (\partial_x z_b)^2} \frac{q}{h} \right) \\ = 2h(\partial_x u)^2 + h \partial_x (g \partial_x(z_b + h) + \partial_x p^{atm}) + h \partial_x \left( \frac{\partial_x z_b}{4 + (\partial_x z_b)^2} (u^2 \partial_{xx}^2 z_b - (g \partial_x(z_b + h) + \partial_x p^{atm}) \partial_x z_b) \right). $$

However, the complexity of the operators in the latter equation made us prefer working with (13).

**Remark 3.** Whether it be from the theoretical or numerical points of view, the equation of interest is

$$ - \nabla_{\text{sgn}} \left( \frac{1}{h} \nabla_{\text{sgn}} \mathbf{Q} \right) = \mathbf{f} \quad (15) $$

for some right hand side f. The well-posedness of Eq. (15) is proven in [4] by means of the Lax-Milgram theorem applied under some smoothness hypotheses on the water height h. In particular, straightforward computations show that

$$ - \int_I \nabla_{\text{sgn}} \left( \frac{1}{h} \nabla_{\text{sgn}} \mathbf{Q} \right) \cdot \tilde{\mathbf{Q}} \, dx = - \left[ \tilde{q} \left( \partial_x(hq) + q_b \partial_x z_b \right) \right]_{\partial I} + \int_I \frac{1}{h} \nabla_{\text{sgn}} \mathbf{Q} \cdot \nabla_{\text{sgn}} \tilde{\mathbf{Q}} \, dx, $$

which explains the present choice of boundary conditions (7a-7b).

Finally, we show that thanks to the correct formulation of diagnostic equations (6e) and (6f), associated to the definition (8b) of $\nabla_{\text{sgn}} \mathbf{X}$ that provides the duality relation (9) we can deduce an associated energy for SGN model. It is crucial from both theoretical and numerical points of view as highlighted in the following section, corresponding to the definition of the numerical method.

**Proposition 4.** Smooth solutions of System (10) satisfy the following energy equality provided $p^{atm}$ and $z_b$ do not depend on time:

$$ \partial_t \left[ h \left( \frac{|\mathbf{X}|^2}{2} + g \left( z_b(x) + \frac{h}{2} \right) + p^{atm}(x) \right) \right] + \partial_x \left[ h u \left( \frac{|\mathbf{X}|^2}{2} + g(z_b(x) + h) + p^{atm}(x) + q \right) \right] = 0. $$

*Proof.* Let us multiply (6b) by u, (6c) by w and (6d) by $\sigma$ so that

$$ \partial_t \left( h \frac{|\mathbf{X}|^2}{2} \right) + \partial_x \left( h u \frac{|\mathbf{X}|^2}{2} \right) + \frac{|\mathbf{X}|^2}{2} \left[ \partial_t h + \partial_x(hu) \right] + \mathbf{X} \cdot \nabla_{\text{sgn}} \mathbf{Q} + h u \partial_x (g(h + z_b) + p^{atm}) = 0. $$

Using the duality relation (9) as well as Eq. (6a) leads to the expected result. □

## 2.2 Splitting strategy at the semi-discrete level

Let us write a semi-discretisation of System (10) based on a classic splitting technique for some time step $\Delta t > 0$, like e.g. in [1, 18, 21, 35, 43]. We first consider the hyperbolic step<fnref n="4" />

$$ \left\{ \begin{array}{l} \displaystyle \frac{h^{n+1/2} - h^n}{\Delta t} + \partial_x(h^n u^n) = 0, \\[2ex] \displaystyle \frac{(h\mathbf{X})^{n+1/2} - (h\mathbf{X})^n}{\Delta t} + \partial_x(h^n u^n \mathbf{X}^n) = \mathbf{S}_{h^n}, \end{array} \right. \quad (16) $$

<footnote n="4">⁴An explicit scheme is presented but an implicit strategy can also be chosen. It does not interfere with the second step.</footnote>
coupled to boundary conditions (7c) and (7d). Any classic numerical method dedicated to the Shallow Water equations (including well-balanced schemes) can be used to solve it.

Then the dispersive step reads

$$ \frac{h^{n+1} - h^{n+1/2}}{\Delta t} = 0, \qquad (17a) $$

$$ \frac{(h\mathbf{X})^{n+1} - (h\mathbf{X})^{n+1/2}}{\Delta t} + \nabla_{\text{sgn}} \mathbf{Q}^{n+1} = 0, \qquad (17b) $$

$$ \nabla_{\text{sgn}} \mathbf{X}^{n+1} = \mathbf{0}, \qquad (17c) $$

which reads as a mixed velocity-pressure Darcy problem. Boundary conditions (7c), (7a) and (7b) are considered. In this paper we investigate a projection-correction approach.

To increase the order of the method, we can consider the incremental method (see e.g. [29]) which consists of the following modified first step

$$ \left\{ \begin{aligned} & \frac{h^{n+1/2} - h^n}{\Delta t} + \partial_x(h^n u^n) = 0, \\ & \frac{(h\mathbf{X})^{n+1/2} - (h\mathbf{X})^n}{\Delta t} + \partial_x(h^n u^n \mathbf{X}^n) + \nabla_{\text{sgn}} \mathbf{Q}^n = \mathbf{S}_{h^n}, \end{aligned} \right. $$

(which requires an initialisation of the pressure unlike the non-incremental version) and the modified second step

$$ \left\{ \begin{aligned} & \frac{h^{n+1} - h^{n+1/2}}{\Delta t} = 0, \\ & \frac{(h\mathbf{X})^{n+1} - (h\mathbf{X})^{n+1/2}}{\Delta t} + \nabla_{\text{sgn}} (\mathbf{Q}^{n+1} - \mathbf{Q}^n) = 0, \\ & \nabla_{\text{sgn}} \mathbf{X}^{n+1} = \mathbf{0}. \end{aligned} \right. $$

This approach does not raise any additional issues, therefore for the sake of simplicity we present the results by considering the splitting strategy (16)-(17).

To make the notation easier, we shall denote from now on: $\mathbf{X} = \mathbf{X}^{n+1}$ and $\mathbf{X}^* = \mathbf{X}^{n+1/2}$. System (17) implies

$$ - \nabla_{\text{sgn}} \left( \frac{1}{h^*} \nabla_{\text{sgn}} \mathbf{Q} \right) = - \frac{\nabla_{\text{sgn}} \mathbf{X}^*}{\Delta t}, \qquad (18) $$

or equivalently, expanding the non-classic operators:

$$ \left\{ \begin{aligned} & 12 \frac{q}{h^*} - h^* \partial_x \left( \frac{\partial_x (h^* q)}{h^*} \right) - 6 \frac{q_b}{h^*} - h^* \partial_x \left( \frac{q_b}{h^*} \partial_x z_b \right) = - \frac{2\sqrt{3}\sigma^* + h^* \partial_x u^*}{\Delta t}, \\ & (4 + (\partial_x z_b)^2) \frac{q_b}{h^*} - 6 \frac{q}{h^*} + \partial_x z_b \frac{\partial_x (h^* q)}{h^*} = - \frac{w^* - u^* \partial_x z_b - \sqrt{3}\sigma^*}{\Delta t}. \end{aligned} \right. \qquad (18') $$

We recover the same operators as in the continuous case – see (13), only the right hand side is modified due to the splitting method. Hence the well-posedness investigated in [4] still holds.

**Remark 4.** The right hand side in (18') is not of order $O(\Delta t^{-1})$ but of order $O(1)$. Indeed, inserting the values from the hyperbolic step (16), we get

$$ \begin{aligned} \frac{2\sqrt{3}\sigma^* + h^* \partial_x u^*}{\Delta t} &= \frac{2\sqrt{3}\sigma^n + h^n \partial_x u^n}{\Delta t} + \mathcal{O}(1) = \mathcal{O}(1), \\ \frac{w^* - u^* \partial_x z_b - \sqrt{3}\sigma^*}{\Delta t} &= \frac{w^n - u^n \partial_x z_b - \sqrt{3}\sigma^n}{\Delta t} + \mathcal{O}(1) = \mathcal{O}(1), \end{aligned} $$

since $2\sqrt{3}\sigma^n + h^n \partial_x u^n = 0$ and $w^n - u^n \partial_x z_b - \sqrt{3}\sigma^n = 0$, provided initial conditions satisfy the divergence constraints, i.e. initial conditions are well-prepared.

Once Equation (18) is solved for $\mathbf{Q}$, the velocity field is updated using (17b)

$$ \mathbf{X} = \mathbf{X}^* - \frac{\Delta t}{h^*} \nabla_{\text{sgn}} \mathbf{Q}. $$
## 2.3 Numerical method at the discrete level

We focus in this section on the discretisation in space of the non-incremental version (16-17). The incremental version does not raise major additional difficulties.

Let us consider a homogeneous Cartesian grid of interval $I = [x_\ell, x_r]$ with mesh size $\Delta x = \frac{x_r-x_\ell}{N}$ for some integer $N > 0$ (see Figure 2).

Standard approaches for Stokes-like problems rely on a staggered grid with velocity fields at the centre of the cells and pressure fields at the interfaces. A similar approach is also followed in [21]. The specific expression of $\nabla_{\text{sgn}} \mathbf{Q}$ in the present problem – see (8b) – induced a different choice, namely a collocated approach with all variables ($\mathbf{X}$ and $\mathbf{Q}$) at the centres of the cells. This choice makes a discrete energy estimate easier to derive.

<figure><img src="image_3.png" /><figcaption>Figure 2: Colocated mesh</figcaption></figure>

### 2.3.1 Hyperbolic step

System (16) consists of the classic shallow water part and two transport equations for $w$ and $\sigma$. A finite-volume numerical scheme reads like this:

$$ \frac{(h\tilde{\mathbf{X}})_i^{n+1/2} - (h\tilde{\mathbf{X}})_i^n}{\Delta t} + \frac{\mathcal{F}_{i+1/2}^n - \mathcal{F}_{i-1/2}^n}{\Delta x} = (0, \mathbf{S}_{h_i^n}^T)^T, $$

where $\tilde{\mathbf{X}} = (1, \mathbf{X}^T)^T$. The eigenvalues of the hyperbolic operator in (16) are $u^n \pm \sqrt{gh^n}$ and $u^n$ (transport equations for $w$ and $\sigma$).

We have considered a high order Polynomial Viscosity Matrix finite-volume method using a second-order MUSCL state reconstruction operator similar to those proposed in [18, 21]. Moreover, a third-order CWENO – [13] reconstruction has been implemented as well. The considered explicit high-order methods are well-balanced for water at rest solutions

$$ \eta = \text{cst}, \quad u = w = \sigma = 0, $$

and linearly $L^\infty$ stable provided the following condition holds

$$ \Delta t \le C_{CFL} \frac{\Delta x}{\max_i (|u_i^n| + \sqrt{gh_i^n})}. \quad (19) $$

Moreover, the schemes are positive preserving for the total water depth in the sense that if $h_i^n > 0$ then $h_i^{n+1} > 0$ for all $i$.

### 2.3.2 Non-hydrostatic step: numerical scheme for the velocity-pressure problem

We could have directly discretised Eqs. (18'). However, in order to recover a discrete energy, we first discretise the velocity-pressure problem (17) and we then deduce a discretisation for (18') mimicking the continuous level. Finite differences are applied in the present work.

Let us first consider the mixed formulation (17). Denote by $\mathcal{H} \in \mathcal{M}_{N,N}(\mathbb{R})$ the diagonal matrix with entries $\mathcal{H}_{i,i} = h_i^*$.
We expect that the discretisation of System (17) has the form

$$ \left( \frac{\bar{\mathcal{H}}/\Delta t}{\mathcal{B}^T} \middle| \frac{\mathcal{B}}{\bar{\bar{\mathbf{\Omega}}}} \right) \left( \frac{\mathbf{X}}{\mathbf{Q}} \right) = \left( \frac{\bar{\mathcal{H}}\mathbf{X}^*/\Delta t - \mathbf{\hat{0}}}{\mathbf{\tilde{\mathbf{0}}}^*} \right). \quad (20) $$

where

* $\bar{\mathcal{H}} \in \mathcal{M}_{3N,3N}(\mathbb{R})$ is block-diagonal with block entries $\mathcal{H} \in \mathcal{M}_{N,N}(\mathbb{R})$;

* $\mathcal{B} \in \mathcal{M}_{3N,2N}(\mathbb{R})$ is an rectangular matrix to be specified;

* $\mathbf{U} = (u_i)_{i \in \{1, \dots, N\}}, \mathbf{W} = (w_i)_{i \in \{1, \dots, N\}}, \mathbf{\Sigma} = (\sigma_i)_{i \in \{1, \dots, N\}}$ and $\mathbf{X} = \begin{pmatrix} \mathbf{U} \\ \mathbf{W} \\ \mathbf{\Sigma} \end{pmatrix}$;

* $\bar{\mathbf{0}}_{ij} = 0$ for $(i, j) \in \{1, \dots, 2N\}^2$;

* $\mathbf{\hat{0}} \in \mathbb{R}^{3N}$ and $\bar{\mathbf{0}} \in \mathbb{R}^{2N}$ account for boundary conditions as detailed below.

The matrix in (20) is symmetric due to the duality relation (9). However, a naive approach where Eqs. (17b) and (17c) are discretised independently from each other would lead to a non-symmetric matrix. Hence we first discretise Eq. (17b) as

$$ \frac{1}{\Delta t} \bar{\mathcal{H}} \mathbf{X} + \mathcal{B} \mathbf{Q} = \frac{1}{\Delta t} \bar{\mathcal{H}} \mathbf{X}^* - \mathbf{\hat{0}}. $$

Once $\mathcal{B}$ is built, $\mathcal{B}^T \mathbf{X} = 0$ should be a consistent discretisation of (17c) up to boundary terms.

More precisely, Equation (17b) is discretised for inner cells:

$$ \frac{h_i^* \mathbf{X}_i}{\Delta t} + \left( \frac{1}{\Delta x} \left( h_{i+1/2}^* \frac{q_{i+1}+q_i}{2} - h_{i-1/2}^* \frac{q_i+q_{i-1}}{2} \right) + q_{bi} (\partial_x z_b)_i \right) = \frac{h_i^* \mathbf{X}_i^*}{\Delta t}, \quad i \in \{2, \dots, N-1\}. $$

which leads to

$$ \mathcal{B} = \begin{pmatrix} B_{11} & B_{12} \\ \bar{\mathbf{0}} & -\mathcal{I}_N \\ -2\sqrt{3}\mathcal{I}_N & \sqrt{3}\mathcal{I}_N \end{pmatrix}. \quad (21) $$

$B_{12}$ is diagonal with entries $(\partial_x z_b)_i$ while $B_{11}$ is tridiagonal:

$$ B_{11} = \frac{1}{\Delta x} \begin{pmatrix} \cdots & \cdots & \cdots & \cdots & \cdots \\ -\frac{h_{i-1/2}^*}{2} & \frac{h_{i+1/2}^*-h_{i-1/2}^*}{2} & \cdots & \cdots & \cdots \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ -\frac{h_{i+1/2}^*}{2} & \frac{h_{i-1/2}^*}{2} & \cdots & \cdots & \cdots \\ \end{pmatrix} \quad (22) $$

Hence, for inner nodes $x_i$:

$$ \begin{aligned} (\mathcal{B}^T \mathbf{X})_i &= - \left[ \frac{1}{2} \left( h_{i-1/2}^* \frac{u_i - u_{i-1}}{\Delta x} + h_{i+1/2}^* \frac{u_{i+1} - u_i}{\Delta x} \right) + 2\sqrt{3}\sigma_i \right], \\ (\mathcal{B}^T \mathbf{X})_{N+i} &= - \left[ w_i - (\partial_x z_b)_i u_i - \sqrt{3}\sigma_i \right], \end{aligned} $$

which is consistent with (17c).
**Boundary conditions.** Let us now specify how boundary conditions are incorporated by means of ghost cells. BC (7a) at $x = x_{1/2} = x_\ell$ is discretised by

$$ \frac{h_1^* q_1 - h_0^* q_0}{\Delta x} = \chi_\ell \implies q_0 = \frac{h_1^* q_1 - \chi_\ell \Delta x}{h_0^*}. $$

Hence, the momentum equation for $u$ is discretised in the first cell by:

$$ \frac{h_1^* u_1}{\Delta t} + \frac{1}{\Delta x} \left[ h_{3/2}^* \frac{q_1 + q_2}{2} - h_{1/2}^* \frac{q_0 + q_1}{2} \right] + q_{b1}(\partial_x z_b)_1 = \frac{h_1^* u_1^*}{\Delta t} \\ \implies \frac{h_1^* u_1}{\Delta t} + \frac{h_{3/2}^*}{2\Delta x} q_2 + \frac{1}{2\Delta x} \left[ h_{3/2}^* - h_{1/2}^* \left( 1 + \frac{h_1^*}{h_0^*} \right) \right] q_1 + q_{b1}(\partial_x z_b)_1 = \frac{h_1^* u_1^*}{\Delta t} - \frac{h_{1/2}^*}{2h_0^*} \chi_\ell. $$

This provides the first line of Matrix $B_{11}$ as well as $\tilde{0}_1 = -\frac{h_{1/2}^*}{2h_0^*} \chi_\ell$. It remains to specify $h_0^*$. Given that no boundary condition is imposed on $h$ at $x = x_\ell$, we set $h_0^* = 2h_1^* - h_2^*$.

We then observe that

$$ (\mathcal{B}^T \mathbf{X})_1 = -\frac{1}{2} \left( h_{3/2}^* \frac{u_2 - u_1}{\Delta x} + h_{1/2}^* \left( 1 + \frac{h_1^*}{h_0^*} \right) \frac{u_1}{\Delta x} \right). $$

We remark that if we discretise BC (7c) by

$$ Q_u = \frac{h_1^* u_1 + h_0^* u_0}{2} \implies \frac{u_1 - u_0}{\Delta x} = \left(1 + \frac{h_1^*}{h_0^*}\right) \frac{u_1}{\Delta x} - \frac{2Q_u}{h_0^* \Delta x}. $$

$$ \text{We deduce } \tilde{0}_1 = \frac{h_{1/2}^*}{h_0^*} \frac{Q_u}{\Delta x}. $$

As for the boundary conditions at $x = x_{N+1/2} = x_r$, BC (7b) is discretised by

$$ \frac{q_N + q_{N+1}}{2} = 0, $$

which means that the momentum equation for $u$ in the last cell becomes

$$ \frac{h_N^* u_N}{\Delta t} - \frac{h_{N-1/2}^*}{2\Delta x} (q_N + q_{N-1}) + q_{bN} (\partial_x z_b)_N = \frac{h_N^* u_N^*}{\Delta t}. $$

This yields the last line of $B_{11}$ and $\tilde{0}_N = 0$. It implies that

$$ (\mathcal{B}^T \mathbf{X})_N = -h_{N-1/2}^* \frac{u_{N-1} - u_N}{2\Delta x} = -\frac{1}{2} \left( h_{N-1/2}^* \frac{u_N - u_{N-1}}{\Delta x} + h_{N+1/2}^* \frac{u_{N+1} - u_N}{\Delta x} \right) $$

provided that $u_{N+1} = u_N$, which is equivalent to imposing

$$ \partial_x u = 0 \text{ at } x = x_r, $$

and $\tilde{0}_N = 0$.

Let us comment the latter result. Imposing $q = 0$ at some boundary is consistent with imposing $\partial_x u = 0$ at the same location. Even if there is no BC required in the hyperbolic problem, it is necessary from the numerical point of view in order to compute the flux at the corresponding interface. We mention that the BC handled in the first step of the algorithm must not be damaged by the second step. Likewise, still for numerical purposes, we impose

$$ \partial_x w = 0, \quad \partial_x \sigma = 0 \quad \text{at } x = x_r. $$
### 2.3.3 Numerical scheme for the pressure problem

For a resolution of the mixed velocity-pressure problem (20) by means of the Uzawa method, see [4]. In the present paper, we rather combine both equations in (20) to obtain

$$ \mathcal{B}^T \overline{\mathcal{H}}^{-1} \mathcal{B} \mathbf{Q} = \mathbf{f} := \frac{\mathcal{B}^T \mathbf{X}^* - \tilde{\mathbf{0}}}{\Delta t} - \mathcal{B}^T \overline{\mathcal{H}}^{-1} \hat{\mathbf{0}}, \quad (23) $$

which is nothing but a consistent discretisation of the velocity-correction approach (18).

**Remark 5.** It is crucial that the initial data satisfy at the discrete level $\mathcal{B}^T \mathbf{X}^* - \tilde{\mathbf{0}} = \mathbf{0}$ in order to prevent the propagation of errors of order $O(\Delta t^{-1})$ as explained in Remark 4

Notice that Matrix $\mathcal{C} = \mathcal{B}^T \overline{\mathcal{H}}^{-1} \mathcal{B}$ has the following structure:

$$ \mathcal{C} = \left( \begin{array}{c|c} C_{11} & C_{12} \\ \hline C_{12}^T & C_{22}^T \end{array} \right) = \left( \begin{array}{c|c} B_{11}^T \mathcal{H}^{-1} B_{11} + 12 \mathcal{H}^{-1} & B_{11}^T \mathcal{H}^{-1} B_{12} - 6 \mathcal{H}^{-1} \\ \hline B_{12}^T \mathcal{H}^{-1} B_{11} - 6 \mathcal{H}^{-1} & B_{12}^T \mathcal{H}^{-1} B_{12} + 4 \mathcal{H}^{-1} \end{array} \right). \quad (24) $$

**Lemma 1.** Matrix $\mathcal{C}$ is symmetric positive-definite.

*Proof.* The symmetry comes from the definition $\mathcal{C} = \mathcal{B}^T \overline{\mathcal{H}}^{-1} \mathcal{B}$. Moreover, we have

$$ \langle \mathcal{C} \mathbf{Q}, \mathbf{Q} \rangle = \langle \mathcal{H}^{-1} \mathcal{B} \mathbf{Q}, \mathcal{B} \mathbf{Q} \rangle \geq 0 $$

since $\mathcal{H}$ is a diagonal matrix with positive entries provided the hyperbolic scheme is positive. Then $\mathcal{B}\mathbf{Q} = \mathbf{0} \Rightarrow \mathbf{Q} = \mathbf{0}$ as $\mathcal{B}$ is one-to-one – see (21).

More precisely, (24) provides the following discrete energy

$$ \langle \mathcal{C} \mathbf{Q}, \mathbf{Q} \rangle = \langle \mathcal{H}^{-1} \mathbf{q}_b, \mathbf{q}_b \rangle + 12 \langle \mathcal{H}^{-1} \left( \mathbf{q} - \frac{\mathbf{q}_b}{2} \right), \mathbf{q} - \frac{\mathbf{q}_b}{2} \rangle + \langle \mathcal{H}^{-1} (B_{11} \mathbf{q} + B_{12} \mathbf{q}_b), B_{11} \mathbf{q} + B_{12} \mathbf{q}_b \rangle, $$

which is consistent with

$$ \int_I \frac{1}{h^*} |\nabla_{\text{sgn}} \mathbf{Q}|^2 dx = \int_I \frac{1}{h^*} \left[ q_b^2 + 12 \left( q - \frac{q_b}{2} \right)^2 + (\partial_x (h^* q) + q_b \partial_x z_b)^2 \right] dx. $$

□

To solve (23), one may apply different strategies. One of them is to apply a direct resolution inverting $\mathcal{C}$. Another variant is to use the Uzawa method. Here, we shall propose to follow an iterative process of Gauss-Seidel type:

1. $C_{11}q^{k+1} = \frac{1}{\Delta t} \mathbf{f} - C_{12}q_b^k$;

2. $C_{22}q_b^{k+1} = \frac{1}{\Delta t} \mathbf{f}_b - C_{12}^T q^{k+1}$.

This method converges due to Lemma 1. This is the latter method which has been chosen in the present work. Matrix $C_{11}$ is more or less the same matrix as for the DAE model (5). The extra additional time then corresponds to the iterative process ($C_{11}$ is factorised once and $C_{22}$ is diagonal).

**Remark 6.** Notice that we could explicitly solve the second block of (23):

$$ q_b = C_{22}^{-1} \left( \frac{1}{\Delta t} \mathbf{f}_b - C_{12}^T \mathbf{q} \right) \Rightarrow (C_{11} - C_{12} C_{22}^{-1} C_{12}^T) \mathbf{q} = \frac{1}{\Delta t} (\mathbf{f} - C_{12} C_{22}^{-1} \mathbf{f}_b), $$

which corresponds to a discretisation of (14).

**Remark 7.** The final numerical scheme proposed here is well-balanced for the water at rest solution, and positive preserving for the total water depth provided the hyperbolic scheme is well-balanced and positive preserving.
<figure><img src="image_4.png" /><figcaption>Figure 3: Multilayer description</figcaption></figure>

# 3 Numerical scheme for the LDNH₂(L) model

## 3.1 Notations and model

Without any assumption on the shallowness of the flow, we can also approximate the Euler equations (1) by means of a "multilayer" model. The flow is split into arbitrary layers of constant height $h_\alpha = \frac{h}{L}$ for some integer $L > 0$ (see Fig. 3). We consider the dispersive layer-averaged model derived in [24] under the following reformulation:

$$ \partial_t h + \partial_x (h \bar{u}) = 0, \quad \bar{u} = \sum_{\alpha=1}^{L} \ell_\alpha \bar{u}_\alpha, $$

and for $\alpha \in \{1, \dots, L\}$

$$ \partial_t(h_\alpha \bar{u}_\alpha) + \partial_x(h_\alpha \bar{u}_\alpha^2 + h_\alpha \bar{q}_\alpha) + \tilde{u}_{\alpha+1/2} \Gamma_{\alpha+1/2} - \partial_x z_{\alpha+1/2} q_{\alpha+1/2} \\ - \tilde{u}_{\alpha-1/2} \Gamma_{\alpha-1/2} + \partial_x z_{\alpha-1/2} q_{\alpha-1/2} = -gh_\alpha \partial_x \eta, $$

$$ \partial_t(h_\alpha \bar{w}_\alpha) + \partial_x(h_\alpha \bar{u}_\alpha \bar{w}_\alpha) + \tilde{w}_{\alpha+1/2} \Gamma_{\alpha+1/2} + q_{\alpha+1/2} - \tilde{w}_{\alpha-1/2} \Gamma_{\alpha-1/2} - q_{\alpha-1/2} = 0, $$

$$ \partial_t(h_\alpha \sigma_\alpha) + \partial_x(h_\alpha \sigma_\alpha \bar{u}_\alpha) = 2\sqrt{3} \left[ \bar{q}_\alpha - \frac{q_{\alpha+1/2} + q_{\alpha-1/2}}{2} \right. \\ \left. - \Gamma_{\alpha+1/2} \left( \frac{h_\alpha \partial_x \bar{u}_\alpha}{12} + \frac{\tilde{w}_{\alpha+1/2} - \bar{w}_\alpha}{2} \right) + \Gamma_{\alpha-1/2} \left( \frac{h_\alpha \partial_x \bar{u}_\alpha}{12} + \frac{\bar{w}_\alpha - \tilde{w}_{\alpha-1/2}}{2} \right) \right], $$

together with diagnostic equations

$$ 2\sqrt{3}\sigma_{\alpha} + h_{\alpha} \partial_{x} \bar{u}_{\alpha} = 0, \quad \alpha \in \{1, \dots, L\}, $$

$$ \bar{w}_{\alpha} - \bar{w}_{\alpha-1} - (\bar{u}_{\alpha} - \bar{u}_{\alpha-1}) \partial_x z_{\alpha-1/2} - \sqrt{3}(\sigma_{\alpha-1} + \sigma_{\alpha}) = 0, \quad \alpha \in \{2, \dots, L\}, $$

$$ w_1 - u_1 \partial_x z_b - \sqrt{3}\sigma_1 = 0, $$
and the boundary condition

$$q_{L+1/2} = 0.$$

The mass transfer term is given by

$$\Gamma_{\alpha+1/2} = \sum_{\beta=\alpha+1}^{L} \partial_x (h_{\beta} (\bar{u}_{\beta} - \bar{\bar{u}})).$$

To close the model, we define

$$\tilde{u}_{\alpha+1/2} = (1 - \gamma_{\alpha+1/2}) \bar{u}_{\alpha+1} + \gamma_{\alpha+1/2} \bar{u}_{\alpha};$$

$$\tilde{w}_{\alpha+1/2} = (1 - \gamma_{\alpha+1/2}) \left( \bar{w}_{\alpha+1} + \frac{h_{\alpha+1}}{2} \partial_x \bar{u}_{\alpha+1} \right) + \gamma_{\alpha+1/2} \left( \bar{w}_{\alpha} - \frac{h_{\alpha}}{2} \partial_x \bar{u}_{\alpha} \right);$$

for any $\gamma_{\alpha+1/2} \in [0, 1]$ such that [24, Prop. 1]

$$\left( \gamma_{\alpha+1/2} - \frac{1}{2} \right) \Gamma_{\alpha+1/2} \geq 0. \qquad (25)$$

Let us recall the properties of this system:

* Let us set $\bar{K}_\alpha = \frac{\bar{u}_\alpha^2 + \bar{w}_\alpha^2 + \sigma_\alpha^2}{2}$. If $(H, \bar{u}_\alpha, \bar{w}_\alpha, \bar{q}_\alpha)$ are smooth solutions to the multilayer model, we have under (25)

$$\partial_t \left( \sum_{\alpha=1}^{L} h_{\alpha} (\bar{K}_{\alpha} + gz_{\alpha} + p^{\text{atm}}(x)) \right) + \partial_x \left( \sum_{\alpha=1}^{L} h_{\alpha} \bar{u}_{\alpha} (\bar{K}_{\alpha} + \bar{q}_{\alpha} + g\eta + p^{\text{atm}}(x)) \right) \le 0. \quad (26)$$

Moreover, if we take $\gamma_{\alpha+1/2} = \frac{1}{2}$, then (26) is an equality.

* The monolayer case $L=1$ reduces to the Serre-Green-Naghdi equations (6).

## 3.2 Differential operators

Let us introduce the operators

$$\mathbf{X}_{\alpha} = \begin{pmatrix} \bar{u}_{\alpha} \\ \bar{w}_{\alpha} \\ \sigma_{\alpha} \end{pmatrix}, \qquad \mathbf{Q}_{\alpha} = \begin{pmatrix} \bar{q}_{\alpha} \\ q_{\alpha-1/2} \end{pmatrix},$$

and

$$\nabla_{ldnh}^{\alpha} \mathbf{Q}_{\alpha} = \begin{pmatrix} \frac{\partial_x (h^* \bar{q}_{\alpha})}{L} - q_{\alpha+1/2} \partial_x z_{\alpha+1/2}^* + q_{\alpha-1/2} \partial_x z_{\alpha-1/2}^* \\ q_{\alpha+1/2} - q_{\alpha-1/2} \\ -2\sqrt{3} \left( \bar{q}_{\alpha} - \frac{q_{\alpha+1/2} + q_{\alpha-1/2}}{2} \right) \end{pmatrix}, \quad \alpha \in \{1, \dots, L\},$$

$$\nabla_{ldnh}^{\alpha} \mathbf{X}_{\alpha} = \begin{pmatrix} 2\sqrt{3}\sigma_{\alpha} + \frac{h_{\alpha}^{*}}{L} \partial_{x} \bar{u}_{\alpha} \\ \bar{w}_{\alpha} - \bar{w}_{\alpha-1} - (\bar{u}_{\alpha} - \bar{u}_{\alpha-1}) \partial_{x} z_{\alpha-1/2}^{*} - \sqrt{3}(\sigma_{\alpha-1} + \sigma_{\alpha}) \end{pmatrix}, \quad \alpha \in \{2, \dots, L\},$$

$$\nabla_{ldnh}^{1} \mathbf{X}_{1} = \begin{pmatrix} 2\sqrt{3}\sigma_{1} + \frac{h_{1}^{*}}{L} \partial_{x} \bar{u}_{1} \\ \bar{w}_{1} - \bar{u}_{1} \partial_{x} z_{b} - \sqrt{3}\sigma_{1} \end{pmatrix}.$$

We get a global duality relation

$$\sum_{\alpha=1}^{L} \mathbf{Q}_{\alpha} \cdot (\nabla_{ldnh}^{\alpha} \mathbf{X}_{\alpha}) = \frac{1}{L} \partial_{x} (h \bar{q}_{\alpha} \bar{u}_{\alpha}) - \sum_{\alpha=1}^{L} \mathbf{X}_{\alpha} \cdot \nabla_{ldnh}^{\alpha} \mathbf{Q}_{\alpha}. \qquad (27)$$
It is worth noticing that

$$ \nabla_{ldnh}^{\alpha} \mathbf{Q}_{\alpha} = \begin{pmatrix} \frac{\partial_x (h^* \bar{q}_{\alpha})}{L} \\ 0 \\ -2\sqrt{3} \bar{q}_{\alpha} \end{pmatrix} + \begin{pmatrix} q_{\alpha-1/2} \partial_x z_{\alpha-1/2}^* \\ -q_{\alpha+1/2} \\ \sqrt{3} q_{\alpha+1/2} \end{pmatrix} - \begin{pmatrix} q_{\alpha+1/2} \partial_x z_{\alpha+1/2}^* \\ -q_{\alpha-1/2} \\ \sqrt{3} q_{\alpha-1/2} \end{pmatrix} + \begin{pmatrix} 0 \\ 0 \\ 2\sqrt{3} q_{\alpha+1/2} \end{pmatrix}, \quad (28a) $$

and

$$ \nabla_{ldnh}^{\alpha} \mathbf{X}_{\alpha} = \left( \bar{w}_{\alpha} - \bar{u}_{\alpha} \partial_x z_{\alpha-1/2}^* - \sqrt{3} \sigma_{\alpha} \right) - \left( \bar{w}_{\alpha-1} - \bar{u}_{\alpha-1} \partial_x z_{\alpha-1/2}^* - \sqrt{3} \sigma_{\alpha-1} \right) - \left( 2\sqrt{3} \sigma_{\alpha-1} \right). \quad (28b) $$

We mention that the divergence operator stated in the present work is not the same as in [24]. We rather provide an equivalent formulation that satisfies the duality relation (27). In particular, $\sigma_\alpha$ is involved instead of $\partial_x \bar{u}_\alpha$ which is important for the order of derivatives.

Let us proceed as in the single layer case with a splitting strategy between hyperbolic terms and non-hydrostatic terms. The latter part reads

$$ \left\{ \begin{aligned} & \frac{h^* \mathbf{X}_\alpha - h^* \mathbf{X}_\alpha^*}{L \Delta t} + \nabla_{ldnh}^\alpha \mathbf{Q}_\alpha = 0, \\ & \nabla_{ldnh}^\alpha \mathbf{X}_\alpha = 0, \end{aligned} \right. \quad (29) $$

which implies

$$ - \nabla_{ldnh}^{\alpha} \left( \frac{\nabla_{ldnh}^{\alpha} \mathbf{Q}_{\alpha}}{h^{*}} \right) = - \frac{1}{L \Delta t} \nabla_{ldnh}^{\alpha} \mathbf{X}^{*}. $$

The latter equation is expanded as, for $\alpha \in \{1, \dots, L\}$

$$ -h^* \partial_x \left( \frac{\partial_x (h^* \bar{q}_\alpha)}{h^*} \right) + 12L^2 \frac{\bar{q}_\alpha}{h^*} = -\frac{L}{\Delta t} \left( 2\sqrt{3}\sigma_\alpha^* + \frac{h^*}{L} \partial_x u_\alpha^* \right) \\ -Lh^* \partial_x \left( \frac{q_{\alpha+1/2} \partial_x z_{\alpha+1/2}^* - q_{\alpha-1/2} \partial_x z_{\alpha-1/2}^*}{h^*} \right) + 6L^2 \frac{q_{\alpha+1/2} + q_{\alpha-1/2}}{h^*}, \quad (30a) $$

then for $\alpha \in \{2, \dots, L\}$

$$ \begin{aligned} & \left(2 - \partial_x z_{\alpha+1/2}^* \partial_x z_{\alpha-1/2}^*\right) \frac{q_{\alpha+1/2}}{h^*} + 2 \left(4 + (\partial_x z_{\alpha-1/2}^*)^2\right) \frac{q_{\alpha-1/2}}{h^*} + \left(2 - \partial_x z_{\alpha-1/2}^* \partial_x z_{\alpha-3/2}^*\right) \frac{q_{\alpha-3/2}}{h^*} \\ & = -\frac{1}{L \Delta t} \left(w_\alpha^* - w_{\alpha-1}^* - (u_\alpha^* - u_{\alpha-1}^*) \partial_x z_{\alpha-1/2}^* - \sqrt{3}(\sigma_\alpha^* + \sigma_\alpha^*)\right) \\ & \quad - \frac{\partial_x z_{\alpha-1/2}^*}{L} \partial_x \left(h^* (\bar{q}_\alpha - \bar{q}_{\alpha-1})\right) + 6 \frac{\bar{q}_{\alpha-1} + \bar{q}_\alpha}{h^*}, \end{aligned} \quad (30b) $$

and finally

$$ \begin{aligned} & (2 - \partial_x z_{3/2}^* \partial_x z_{1/2}^*) \frac{q_{3/2}}{h^*} + (4 + (\partial_x z_{1/2}^*)^2) \frac{q_{1/2}}{h^*} \\ & = -\frac{1}{L \Delta t} \left( w_1^* - u_1^* \partial_x z_b - \sqrt{3} \sigma_1^* \right) - \frac{\partial_x z_b}{L} \frac{\partial_x (h^* \bar{q}_1)}{h^*} + 6 \frac{\bar{q}_1}{h^*}. \end{aligned} \quad (30c) $$

Let us remark is that the differential operator for $\bar{q}_\alpha$ in (30a) is independent from index $\alpha$. Moreover, up to a coefficient $L^2$, this operator is the same as in the monolayer case – see (18').

### 3.3 Structure of the fully discretised equations

Our goal is to produce an algorithm that only relies on tools designed for the monolayer case. Let us split the gradient matrix $\mathcal{B} = \mathcal{B}^{1/2}$ from (21) into:

$$ \mathcal{B}^{\alpha-1/2} = \mathcal{B}_1 + \mathcal{B}_2^{\alpha-1/2}, \text{ where } \mathcal{B}_1 = \begin{pmatrix} B_{11} & \bar{0} \\ \bar{0} & \bar{0} \\ -2\sqrt{3} \mathcal{I}_N & \bar{0} \end{pmatrix} \quad \text{and} \quad \mathcal{B}_2^{\alpha-1/2} = \begin{pmatrix} \bar{0} & B_{12}^{\alpha-1/2} \\ \bar{0} & -\mathcal{I}_N \\ \bar{0} & \sqrt{3} \mathcal{I}_N \end{pmatrix}, $$
• $C_{1,1} = C^{1/2}$;

• For $\alpha \in \{2, \dots, L\}$:

$$C_{\alpha,\alpha} = C^{\alpha-1/2} + (\mathcal{R} - \mathcal{B}_2^{\alpha-1/2})^T \mathcal{H}^{-1} (\mathcal{R} - \mathcal{B}_2^{\alpha-1/2})^T = \left( \frac{C_{11}}{(C_{12}^{\alpha-1/2})^T} \middle| \frac{C_{12}^{\alpha-1/2}}{2C_{22}^{\alpha-1/2}} \right)$$

and

$$C_{\alpha-1,\alpha} = C_{\alpha,\alpha-1}^{T} = (\mathcal{B}^{\alpha-3/2})^{T} \mathcal{H}^{-1} (\mathcal{R} - \mathcal{B}_{2}^{\alpha-1/2}) = \left( \frac{\bar{\mathbf{o}}}{\bar{\mathbf{o}}} \middle| \frac{-6\mathcal{H}^{-1} - B_{11}^{T}\mathcal{H}^{-1}B_{12}^{\alpha-1/2}}{2\mathcal{H}^{-1} - (B_{12}^{\alpha-3/2})^{T}\mathcal{H}^{-1}B_{12}^{\alpha-1/2}} \right).$$

(31) is a consistent discretisation of (30).

## 3.4 Iterative scheme

To solve the $2NL \times 2NL$ linear system (31), we apply an alternating direction-type method. More precisely, it amounts to solving iteratively:

**x-direction** For each layer $\alpha$, (30a) in $\bar{q}_\alpha$ knowing $q_{\alpha+1/2}$ which corresponds to a linear system with matrix $C_{11}$ given in (24) (the same for all layers which requires a single factorisation for all iterates at each time step):

$$C_{11} \mathbf{q}_{\alpha}^{p+1} = \frac{1}{L \Delta t} \left( B_{11}^{T} \mathbf{U}_{\alpha}^{*} - 2\sqrt{3} \mathbf{\Sigma}_{\alpha}^{*} \right) - C_{12}^{\alpha-1/2} \mathbf{q}_{\alpha-1/2}^{p};$$

**z-direction** For each node $x_i$, (30b-30c) in $q_{\alpha+1/2}$ knowing $\bar{q}_\alpha$. It is a tridiagonal system with matrix $S^{(i)}$ solved by means of the Thomas' algorithm [48]. Matrix $S^{(i)}$ is a $L \times L$ tridiagonal symmetric positive-definite matrix with

• $S_{1,1}^{(i)} = 4 + (\partial_x z_{1/2}^*)^2_i$;

• For $\alpha \ge 2$, $S_{\alpha,\alpha}^{(i)} = 2(4 + (\partial_x z_{\alpha-1/2}^*)^2_i)$ and $S_{\alpha,\alpha-1}^{(i)} = S_{\alpha-1,\alpha}^{(i)} = 2 - (\partial_x z_{\alpha-1/2}^*)_i (\partial_x z_{\alpha-3/2}^*)_i$.

Indeed, we check that

$$\langle S^{(i)} \mathbf{x}, \mathbf{x} \rangle = 2x_1^2 + (2 + (\partial_x z_{L-1/2})_i^2) x_L^2 + \sum_{\alpha=2}^{L-1} \left[ 4x_\alpha^2 + 2(x_\alpha + x_{\alpha-1})^2 + ((\partial_x z_{\alpha-1/2})_i^* x_\alpha - (\partial_x z_{\alpha-3/2})_i^* x_{\alpha-1})^2 \right] \ge 0.$$

This strategy is equivalent to solving the following system by means of a Gauss-Seidel iterative procedure:

$$\left( \begin{array}{c|c} C_{11} & \bar{\mathbf{o}} \\ \vdots & \vdots \\ \bar{\mathbf{o}} & C_{11} \\ \hline E^T & \end{array} \right) \left( \begin{array}{c} E \\ \mathcal{S}^{(1)} \\ \vdots \\ \mathcal{S}^{(N)} \end{array} \right)$$

where unknowns have been re-labelled as

$$\left( \begin{array}{c} \hat{q}_{\alpha=1} \\ \vdots \\ \hat{q}_{\alpha=L} \\ \hline \check{q}_{i=1} \\ \vdots \\ \check{q}_{i=N} \end{array} \right), \text{ with } \hat{q}_{\alpha} = \left( \begin{array}{c} q_{\alpha,1} \\ \vdots \\ q_{\alpha,N} \end{array} \right) \text{ and } \check{q}_i = \left( \begin{array}{c} q_{1/2,i} \\ \vdots \\ q_{L-1/2,i} \end{array} \right).$$
**Remark 8.** It must be underlined that this algorithm can be easily parallelised insofar as each direction (x and z) involves a blockwise diagonal matrix.

# 4 Numerical simulations

Let us recall that the stability condition is prescribed by the hyperbolic part of the splitting strategy – see (19) with a CFL number $C_{CFL}$ to be specified.

The numerical schemes presented in this paper, namely the resolution of (23) for $NHML_2(L = 1)$ and the resolution of (31) for $NHML_2(L), L \ge 1$, are assessed by means of 4 classic test cases and compared to numerical results obtained with the DAE model (5) and its multilayer counterpart described in [24]:

* Propagation of the Euler solitary wave;

* Comparisons with 3 sets of experimental data.

Let us mention that in order to ensure the robustness of the algorithm for wet/dry transitions, the term $\frac{1}{h_i^n}$ is approximated by [38]

$$ \frac{h_i^n \sqrt{2}}{\sqrt{(h_i^n)^4 + \max\{(h_i^n)^4, \varepsilon\}}}, $$

for some threshold $\varepsilon > 0$, as it is usually done in shallow waters (see [32]).

## 4.1 Convergence test

Let us start by a convergence test for both models, $LDNH_0$ and $LDNH_2$, in order to assess the numerical strategy and the code used. To do so, we consider the propagation of a soliton (specific to each model) in a rectangular channel with constant topography. This solitary wave is given by

$$ \begin{aligned} h(t, x) &= H_* + A \operatorname{sech}^2 \left( \frac{1}{H_*} \sqrt{\frac{A\gamma}{2(A + H_*)}} (x - ct) \right), \quad u(t, x) = c \left( 1 - \frac{H_*}{h(t, x)} \right), \\ w(t, x) &= -\frac{cA}{h(t, x)} \sqrt{\frac{A\gamma}{2(A + H_*)}} \operatorname{sech}^3 \left( \frac{1}{H_*} \sqrt{\frac{A\gamma}{2(A + H_*)}} (x - ct) \right) \sinh \left( \frac{1}{H_*} \sqrt{\frac{A\gamma}{2(A + H_*)}} (x - ct) \right) \end{aligned} \quad (32) $$

where $A$ and $H_*$ are constant fixed values, $c = \sqrt{g(A + H_*)}$ and

* $\gamma = 2$ in the case of $LDNH_0$,

* $\gamma = 3/2$ in the case of $LDNH_2$.

Here we set $H_* = 1$, $A = 0.1$, $g = 1$ and $z_b = 0$.

The propagation of a solitary wave over a long distance is a standard assessment of stability and conservative properties of numerical schemes for Boussinesq-type equations [45, 46, 49]. A solitary wave propagates at constant speed and without change of shape over a horizontal bottom.

The domain is $[-20, 20]$. We perform the simulation with different numbers of volume cells at time $t = 0.2$ with a second-order scheme. The results are compared to the reference solution and the errors are shown in Tables 1 and 2. These results show the convergence towards the reference solution at second-order.
<table><caption>Table 1: Test 4.1 – L¹ errors with the LDNH₀-soliton and numerical orders of accuracy.</caption><thead><tr><th rowspan="2">No. of cells</th><th colspan="2">h</th><th colspan="2">hu</th><th colspan="2">hw</th></tr><tr><th>L<sup>1</sup> error</th><th>order</th><th>L<sup>1</sup> error</th><th>order</th><th>L<sup>1</sup> error</th><th>order</th></tr></thead><tbody><tr><td>50</td><td>3.04e-03</td><td>0.00</td><td>1.81e-03</td><td>0.00</td><td>2.17e-03</td><td>0.00</td></tr><tr><td>100</td><td>9.26e-04</td><td>1.72</td><td>5.17e-04</td><td>1.81</td><td>5.97e-04</td><td>1.86</td></tr><tr><td>200</td><td>2.57e-04</td><td>1.85</td><td>1.39e-04</td><td>1.89</td><td>1.88e-04</td><td>1.67</td></tr><tr><td>400</td><td>6.96e-05</td><td>1.88</td><td>4.22e-05</td><td>1.72</td><td>8.44e-05</td><td>1.15</td></tr></tbody></table>

<table><caption>Table 2: Test 4.1 – L¹ errors with the LDNH₂-soliton and numerical orders of accuracy.</caption><thead><tr><th rowspan="2">No. of cells</th><th colspan="2">h</th><th colspan="2">hu</th><th colspan="2">hw</th></tr><tr><th>L<sup>1</sup> error</th><th>order</th><th>L<sup>1</sup> error</th><th>order</th><th>L<sup>1</sup> error</th><th>order</th></tr></thead><tbody><tr><td>50</td><td>1.21e-02</td><td>0.00</td><td>1.72e-02</td><td>0.00</td><td>1.22e-02</td><td>0.00</td></tr><tr><td>100</td><td>4.51e-03</td><td>1.42</td><td>4.02e-03</td><td>2.10</td><td>3.28e-03</td><td>1.89</td></tr><tr><td>200</td><td>1.33e-03</td><td>1.76</td><td>1.12e-03</td><td>1.85</td><td>8.50e-04</td><td>1.95</td></tr><tr><td>400</td><td>3.80e-04</td><td>1.80</td><td>3.24e-04</td><td>1.79</td><td>2.39e-04</td><td>1.83</td></tr></tbody></table>

## 4.2 Solitary wave propagation over reefs

A test case propagating a solitary wave over an idealised fringing reef assesses the ability of the model to handle nonlinear dispersive waves, breaking waves and bore propagation. The test configuration includes a fore reef, a flat reef, and an optional reef crest to represent fringing reefs commonly found in a tropical environment. Figure 4 shows a sketch of the laboratory experiments carried out at the O.H. Hinsdale Wave Research Laboratory of Oregon State University. See for instance [50] for more details. The 1D domain [0, 45] is discretised with $\Delta x = 0.045$ m.

<figure><img src="image_5.png" /><figcaption>Figure 4: Test 4.2 – sketch of the topography</figcaption></figure>

A solitary wave of amplitude 0.5 m is placed at point $x = 10$ m. Finally $CFL = 0.9$ and $g = 9.81 \text{ m} \cdot \text{s}^{-2}$. Free outflow boundary conditions are imposed.

Figure 5 shows snapshots at different times, $t\sqrt{g/H_*} = t_0$ where $H_* = 1$ m. Comparisons between experimental and simulated data allow to validate the numerical approach presented in this paper. Results are shown for $LDNH_0(L = 1)$ and $LDNH_2(L = 1)$ models. The water rushes over the flat reef without producing a pronounced bore-shape. The simulation also captures the offshore component of the rarefaction falls, exposing the reef edge, below the initial water level. The simulations match with experimental data and the $LDNH_2$ provides slightly better results.

## 4.3 Wave propagation over a submerged bar

The Dingemans experiment [15] of plunging breaking periodic waves over a submerged bar is considered. This case allows to study frequency dispersion characteristics and non-linear interactions. As waves propagate over a submerged bar, multiple phenomena occur, like the appearance of higher harmonics.
<figure><img src="image_6.png" /><figcaption>Figure 5: Test 4.2 – comparison between experimental data (red points) and numerical results (solid blue line for $LDNH_0$ and dashed pink line for $LDNH_2$) at times $t\sqrt{g/H_*} = 0, 7.8, 10, 13, 17, 20.5$ s</figcaption></figure>

The 1D domain $[0, 30]$ is discretised with $\Delta x = 0.005$ and the bathymetry is defined on Fig. 6. Locations of the measurement points are specified in Table 3. The CFL number is set to $C_{CFL} = 0.9$ and the gravity field to $g = 9.81 \text{ m} \cdot \text{s}^{-2}$. We run the numerical test from the “lake at rest” steady state as an initial condition. Boundary conditions correspond to free outflow at $x = 30$ and a sinusoidal wave train for $\eta$ generated at $x = 0$. This is done as in [18] imposing in a relaxation zone:

$$ \eta_{\ell}(t) = A \sin \left( \frac{2\pi}{T} t \right), $$

where $A = 0.01$ and $T = 2.02$ denote resp. amplitude and period.

This test produces, up to the front slope, waves with wavenumbers $k \approx 0.63/H_*$ and $k \approx 1.58H_*$ respectively, where $H_* = 0.4$ is the typical depth. Fig. 7 shows numerical results of time series of the free surface for Model $LDNH_0(L)$, for $L \in \{1, 2, 4\}$, while Fig. 8 concerns $LDNH_2(L)$ and Fig. 9 shows comparisons between $LDNH_0(L = 4)$ and $LDNH_2(L = 4)$.

<figure><img src="image_7.png" /><figcaption>Figure 6: Test 4.3 – sketch of the bathymetry</figcaption></figure>
<table><caption>Table 3: Test 4.3 – wave gauge locations</caption><thead><tr><th>Label</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><td>Location</td><td>10.5</td><td>12.5</td><td>13.5</td><td>14.5</td><td>15.7</td><td>17.3</td><td>19</td><td>21</td></tr></tbody></table>

<table><caption>Table 4: Test 4.4 – position of the wave probes.</caption><thead><tr><th>A</th><th>x<sub>1</sub></th><th>x<sub>2</sub></th><th>x<sub>3</sub></th><th>x<sub>4</sub></th><th>x<sub>5</sub></th></tr></thead><tbody><tr><td>0.289</td><td>23.520</td><td>23.735</td><td>23.990</td><td>24.210</td><td>24.448</td></tr></tbody></table>

Good agreements with experimental data are observed for all models up to Gauge #4. Beyond the bar, higher harmonics are released which explains discrepancies. We recover observations from the literature, such as [39] where σ-coordinates are used, or [18] where an enhanced two-layer version of the non-hydrostatic pressure system $LDNH_0$ is used. The results in [11] with a three-parameter Green-Naghdi model optimised for uneven bottoms, show the same level of agreement. Here, we would like to stress the ability of the proposed models to deal with a wide range of dispersive waves. The main difference between $LDNH_0(4)$ and $LDNH_2(4)$, as pictured on Fig. 9, can be seen for gauges #6 to #8 where the $LDNH_2$ model is more accurate.

## 4.4 Shoaling of a solitary wave on a plane beach

We finally consider the shoaling of a solitary wave on a beach with a constant slope (1 : 30) as described by Guibourg in [30] and then investigated in [6, 16].

A sketch of the geometry is described on Fig. 10. The initial condition is a solitary wave at location $x = 10$ with amplitude $A = 0.298$, as described in [12, 17]. The computational domain $\Omega = [0, 27.5]$ is divided into cells of length $\Delta x = 0.01$. Free-outflow boundary conditions are considered and the CFL number is set to $C_{CFL} = 0.9$.

We compare Models $LDNH_0(L)$ and $LDNH_2(L)$ for $L \in \{1, 2, 4\}$. Some temporal series of the free-surface elevation are measured at various locations (see Table 4) and compared with the corresponding numerical results.

Results are shown on Figs. 11. Numerical outputs for $A = 0.289$ predict the shoaling phenomenon during the wave run-up satisfactorily. The output clearly shows better performance for Model $LDNH_2$, which is expected according to the linear dispersion relation of the continuum models.

## 5 Conclusion

This paper deals with the numerical approach for the models $LDNH_0$ and $LDNH_2$ introduced in [24]. The main objective is to compare their accuracy when applied to different standard test case scenarios. $LDNH_2$ model presented in [24] may be seen as a multilayer extension of the Serre – Green-Naghdi equations. The model was derived from Euler equations assuming linear and quadratic vertical profiles for the vertical velocity and pressure. $LDNH_0$ model presented in [24] differs from $LDNH_2$ on the assumption of a linear vertical profile for pressure. One of the most attractive properties of the models relies on the increasing accuracy of the linear dispersion relation as the number of layers increases, whereas the model $LDNH_2$ shows better accuracy than $LDNH_0$ for the same number of layers.

The $LDNH_0$ model is solved by using a projection technique similar to the one introduced in [18]. The extension of this technique to $LDNH_2$ is not straightforward. We have proposed a numerical method based on this projection technique to approximate the solution of the $LDNH_2$ model.

The complexity of the model $LDNH_2$ requires the design of an efficient strategy to solve it numerically for
<figure><img src="image_8.png" /><figcaption>Figure 7: Test 4.3 – comparison of experiment data (red points) and simulated ones with the model LDNH₀ setting 1 layer (green), 2 layers (blue) and 4 layers (black)</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Test 4.3 – comparison of experiment data (red points) and simulated ones with the model $LDNH_2$ setting 1 layer (green), 2 layers (blue) and 4 layers (black).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Test 4.3 – comparison of experiment data (red points) and simulated ones with the model $LDNH_0$ (blue) and $LDNH_2$ (black) setting 4 layers.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Test 4.4 – sketch of the bathymetry used for the shoaling of a solitary wave test problem.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Test 4.4 with A = 0.289 – blue lines indicate the numerical solution with models $LDNH_0$ (a) and $LDNH_2$ (b) for different layers; red points represent the experimental measurements at the probe positions</figcaption></figure>
an increasing number of layers. To this aim, we have exploited here a duality relation at the continuous
level. This allows to design an algorithm relying on an iterative process that solves a monolayer case in
each iteration.

In particular, the algorithm may be decomposed into two different problems: The first one corresponds to a discrete parabolic problem per layer and the second to a tridiagonal linear system at each point of the horizontal discretization that couples each layer. Moreover, (i) the matrix of the linear system of the parabolic problem per layer is the same for all layers, and (ii) each tridiagonal linear system at each point of the horizontal discretization is independent. Then, we can observe that (i) since we have the same matrix for all layers, we may consider, for example, an LU factorization to diminish the computational time in all layers and (ii) this implies that, although in the paper this is not done, the proposed technique is easily parallelizable, especially when the number of layers increases.

Moreover, the final numerical scheme proposed here is also high-order, well-balanced for the water at
rest solution, and positive preserving for the total water depth. That results in an efficient and robust
numerical scheme, even in the presence of wet/dry transitions. It is worth mentioning that, while an
extensive literature is dedicated to the numerical resolution of shallow water flows such as the Depth-
Averaged Euler equations or the Serre – Green-Naghdi equations, this is the first attempt, up to our
knowledge, of designing a robust and efficient numerical strategy for the multilayer extension of the SGN
equations.

The proposed numerical scheme has been carefully validated, showing the second order of accuracy and comparing it with available experimental data. The obtained results exhibit an excellent fit with the experiments and show that the proposed strategy is well-suited for most coastal processes: wave propagation, shoaling of the waves, run-up of waves onto a beach, higher dispersive harmonic waves, among others.

In the numerical tests presented in the paper, we can also see that for a fixed number of layers, the $LDNH_2$ is in many situations more accurate than $LDNH_0$. This result shows that the $LDNH_2$ is of interest since the accuracy of both models has been compared by using a generalization of classic projection techniques. One should take into account that $LDNH_2$ has twice more pressure unknowns with respect to $LDNH_0$. Therefore, if we apply a projection method to approximate the solution of $LDNH_2$, it results in a more expensive algorithm, from the computational point of view, than for $LDNH_0$. Hence, a fair comparison of both models from the computational point of view is beyond the scope of this paper.

Nevertheless, it is worth mentioning that this paper was aimed to analyse and propose a projection
technique as a starting point for designing efficient and robust numerical methods of high-order for mul-
tilayer non-hydrostatic systems. As further research, it will be interesting to investigate the development
of other numerical strategies for the $LDNH_2$ model and compare it with the one proposed here. For
example, it may be of interest to consider relaxation techniques or rapid numerical methods and pseudo-
compressibility approximations, as the ones proposed for one-layer dispersive shallow models in [19], [23]
and [5], respectively.

# 6 Acknowledgements

This research has been partially supported by the Spanish Government and FEDER through the coordi-
nated Research project RTI2018-096064-B-C1 and RTI2018-096064-B-C2, and the Andalusian Government
Research project UMA18-FEDERJA-161.
References

[1] N. Aïssiouene, M.-O. Bristeau, E. Godlewski, A. Mangeney, and C. Parés. A two-dimensional method for a family of dispersive shallow water models. *The SMAI journal of computational mathematics*, 6:187-226, 2020.

[2] N. Aïssiouene, M.-O. Bristeau, E. Godlewski, and J. Sainte-Marie. A combined finite volume-finite element scheme for a dispersive shallow water system. *Netw. Heterog. Media*, 11(1):1–27, 2016.

[3] S. Allgeyer, M.-O. Bristeau, D. Froger, R. Hamouda, V. Jauzein, A. Mangeney, J. Sainte-Marie, F. Souillé, and M. Vallée. Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface. *ESAIM Math. Model. Numer. Anal.*, 53(6):1981–2024, 2019.

[4] K. Benyo, A. Charhabil, M.-A. Debyaoui, and Y. Penel. Simulation of complex free surface flows. Proc. of the 2019 session of CEMRACS.

[5] A. Bonnet-Ben Dhia, M. Bristeau, E. Godlewski, S. Imperiale, A. Mangeney, and J. Sainte-Marie. Pseudo-compressibility, dispersive model and acoustic waves in shallow water flows. *Preprint. HAL Id: hal-02493518*, 2020.

[6] P. Bonneton, F. Chazel, D. Lannes, F. Marche, and M. Tissier. A splitting approach for the fully nonlinear and weakly dispersive Green–Naghdi model. *J. Comput. Phys.*, 230(4):1479–1498, 2011.

[7] F. Bouchut. *Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources*. Birkhauser, 2004.

[8] J. Boussinesq. Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond. *J. Math. Pures Appl.*, pages 55-108, 1872.

[9] M.-O. Bristeau, B. Di Martino, C. Guichard, and J. Sainte-Marie. Layer-averaged Euler and Navier-Stokes equations. *Commun. Math. Sci.*, 15(5):1221–1246, 2017.

[10] M.-O. Bristeau and J. Sainte-Marie. Derivation of a non-hydrostatic shallow water model; Comparison with Saint-Venant and Boussinesq systems. *Discrete Contin. Dyn. Syst. Ser. B*, 10(4), 2008.

[11] F. Chazel, D. Lannes, and F. Marche. Numerical simulation of strongly nonlinear and dispersive waves using a Green–Naghdi model. *J. Sci. Comput.*, 48(1):105–116, 2011.

[12] D. Clamond and D. Dutykh. https://es.mathworks.com/matlabcentral/fileexchange/39189-solitary-water-wave, 2012.

[13] I. Cravero, G. Puppo, M. Semplice, and G. Visconti. CWENO: uniformly accurate reconstructions for balance laws. *Math. Comp.*, 87(312):1689–1719, 2018.

[14] A. Barré 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.*, 73:237-240, 1871.

[15] M.W. Dingemans. Comparison of computations with boussinesq-like models and laboratory measurements. *Report H-1684.12, 32, Delft Hydraulics*, 1994.

[16] A. Duran and F. Marche. Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations. *Commun. Comput. Phys.*, 17(3):721-760, 2015.

[17] D. Dutykh and D. Clamond. Efficient computation of steady solitary gravity waves. *Wave Motion*, 51(1):86–99, 2014.

[18] C. Escalante, M.J. Castro, E.D. Fernández-Nieto, and T. Morales de Luna. An efficient two-layer non-hydrostatic approach for dispersive water waves. *J. Sci. Comput.*, 79:273–320, 2019.
[19] C. Escalante, M. Dumbser, and M. J. Castro. An efficient hyperbolic relaxation system for dispersive non-hydrostatic water waves and its solution with high order discontinuous Galerkin schemes. *J. Comput. Phys.*, 394:385-416, 2019.

[20] C. Escalante and T. Morales de Luna. A general non-hydrostatic hyperbolic formulation for boussinesq dispersive shallow flows and its numerical approximation. *Journal of Scientific Computing*, 83(3):62, June 2020.

[21] C. Escalante, T. Morales de Luna, and M. J. Castro. Non-hydrostatic pressure shallow flows: GPU implementation using finite volume and finite difference scheme. *Applied Mathematics and Computation*, 338:631-659, December 2018.

[22] N. Favrie and S. Gavrilyuk. A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves. *Nonlinearity*, 30(7), 2017.

[23] N. Favrie and S. Gavrilyuk. A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves. *Non-linearity*, 30(7):2718-2736, 2017.

[24] 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. *Commun. Math. Sci.*, 16(05):1169-1202, 2018.

[25] S. Ferrari and F. Saleri. A new two-dimensional Shallow Water model including pressure effects and slow varying bottom topography. *ESAIM: Math. Model. Numer. Anal.*, 38(2):211-234, 2004.

[26] J.-F. Gerbeau and B. Perthame. Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation. *Discrete Contin. Dyn. Syst. Ser. B*, 1(1):89-102, 2001.

[27] E. Godlewski and P.-A. Raviart. *Numerical approximation of hyperbolic systems of conservation laws*. Number 118 in Applied Mathematical Sciences. Springer-Verlag, 1996.

[28] A.E. Green and P.M. Naghdi. A derivation of equations for wave propagation in water of variable depth. *J. Fluid Mech.*, 78(02):237-246, 1976.

[29] J.-L. Guermond, P. Minev, and J. Shen. An overview of projection methods for incompressible flows. *Comput. Methods App. Mech. Eng.*, 195(44-47):6011-6045, 2006.

[30] S. Guibourg. *Modélisation numérique et expérimentale des houles bidimensionnelles en zone cotière*. PhD thesis, Université Joseph Fourier, 1994.

[31] J.T. Kirby. Boussinesq models and their application to coastal processes across a wide range of scales. *J. Waterway, Port, Coastal, Ocean Eng.*, 142(6):03116005, 2016.

[32] Alexander Kurganov and Guergana Petrova. A second-order well-balanced positivity preserving central-upwind scheme for the saint-venant system. *Commun. Math. Sci.*, 5(1):133-160, 03 2007.

[33] D. Lannes. *The Water Waves Problem: Mathematical Analysis and Asymptotics*, volume 188. Amer. Math. Soc., 2013.

[34] D. Lannes and P. Bonneton. Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation. *Phys. Fluids*, 21(1):016601, 2009.

[35] D. Lannes and F. Marche. A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations. *J. Comput. Phys.*, 282:238-268, 2015.

[36] O. Le Métayer, S. Gavrilyuk, and S. Hank. A numerical scheme for the Green-Naghdi model. *J. Comput. Phys.*, 229:2034-2045, 2010.

[37] R.J. LeVeque. *Finite volume methods for hyperbolic problems*. Number 31. Cambridge Univ. Press, 2002.
[38] 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. *J. Comput. Phys.*, 374:213–236, 2018.

[39] G. Ma, F. Shi, and J.T. Kirby. Shock-capturing non-hydrostatic model for fully dispersive surface wave processes. *Ocean Model.*, 43:22–35, 2012.

[40] P.A. Madsen, R. Murray, and O.R. Sørensen. A new form of the Boussinesq equations with improved linear dispersion characteristics. *Coastal Eng.*, 15(4):371–388, 1991.

[41] F. Marche. Derivation of a new two-dimensional viscous shallow water model with varying topography, bottom friction and capillary effects. *Eur. J. Mech. B Fluids*, 26(1):49–63, 2007.

[42] O. Nwogu. Alternative form of Boussinesq equations for nearshore wave propagation. *J. Waterway, Port, Coastal, Ocean Eng.*, 119(6):618–638, 1993.

[43] M. Parisot. Entropy-satisfying scheme for a hierarchy of dispersive reduced models of free surface flow. *Int. J. Numer. Methods Fluids*, 91(10):509–531, 2019.

[44] D.H. Peregrine. Long waves on a beach. *J. Fluid Mech.*, 27(04):815–827, 1967.

[45] M. Ricchiuto and A.G. Filippini. Upwind residual discretization of enhanced boussinesq equations for wave propagation over complex bathymetries. *J. Comput. Phys.*, 271:306–341, 2014.

[46] V. Roeber, K.F. Cheung, and M.H. Kobayashi. Shock-capturing boussinesq-type model for nearshore wave processes. *Coastal Eng.*, 57:407–423, 2010.

[47] F. Serre. Contribution à l'étude des écoulements permanents et variables dans les canaux. *La Houille Blanche*, (6):830–872, 1953.

[48] L. Thomas. Elliptic problems in linear differential equations over a network. Technical report, Watson scientific computing laboratory, Columbia Univ., NY, 1949.

[49] Y. Yamazaki, Z. Kowalik, and K.F. Cheung. Depth-integrated, non-hydrostatic model for wave breaking and run-up. *Int. J. Numer. Methods Fluids*, 61(5):473–497, 2009.

[50] Yu Yao, Tiancheng He, Zhengzhi Deng, Long Chen, and Huiqun Guo. Large eddy simulation modeling of tsunami-like solitary wave processes over fringing reefs. *Natural Hazards and Earth System Sciences*, 19:1281–1295, 06 2019.