# A geometrical Green-Naghdi type system for dispersive-like waves in prismatic channels

Sergey Gavrilyuk, Mario Ricchiuto

► To cite this version:

Sergey Gavrilyuk, Mario Ricchiuto. A geometrical Green-Naghdi type system for dispersive-like waves in prismatic channels. Journal of Fluid Mechanics, 2025, 1014, pp.A5. 10.1017/jfm.2025.10274. hal-04889204

HAL Id: hal-04889204

https://inria.hal.science/hal-04889204v1

Submitted on 15 Jan 2025

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# A geometrical Green-Naghdi type system for dispersive-like waves in prismatic channels

Sergey Gavrilyuk\* and Mario Ricchiuto†

December 17, 2024

## ABSTRACT

We consider 2D free surface gravity waves in prismatic channels with bathymetric variations uniquely in the transverse direction. Starting from the Saint-Venant equations (shallow water equations) we derive a 1D transverse averaged model describing dispersive effects solely related to variations of the channel topography. These effects have been demonstrated in Chassagne et al. JFM 2019 to be predominant in the propagation of bores with Froude numbers below a critical value of about 1.15. The model proposed is fully nonlinear, Galilean invariant, and admits a variational formulation under natural assumptions about the channel geometry. It is endowed with an exact energy conservation law, and admits exact travelling wave solutions. Our model generalizes and improves the linear equations proposed by Chassagne et al. JFM 2019, as well as Quezada de Luna and Ketcheson JFM 2021. The system is recast in two useful forms appropriate for its numerical approximations, whose properties are discussed. Numerical results allow to verify against analytical solutions the implementation of these formulations, and validate our model against fully 2D nonlinear shallow water simulations, as well as the famous experiments by Treske J. Hyd. Res. 1994.

## Contents

<table><tr><td><b>1 Introduction</b></td><td><b>2</b></td></tr><tr><td><b>2 Equations averaged in y-direction</b></td><td><b>3</b></td></tr><tr><td>2.1 Generalities and asymptotic expansion</td><td>4</td></tr><tr><td>2.2 Small-scale geometrical effects on the section average flow</td><td>5</td></tr><tr><td>2.3 Expressions for the geometric parameter for the triangular and trapezoidal sections</td><td>8</td></tr><tr><td><b>3 Properties of the geometrical Green-Naghdi type system</b></td><td><b>8</b></td></tr><tr><td>3.1 Galilean invariance</td><td>8</td></tr><tr><td>3.2 Lagrangian structure and energy equation</td><td>9</td></tr><tr><td>3.3 Linear dispersion properties</td><td>9</td></tr><tr><td>3.4 Travelling wave solutions: solitary waves</td><td>10</td></tr><tr><td>3.5 Travelling waves solutions: a composite solution</td><td>12</td></tr><tr><td><b>4 System reformulations and numerical approximation</b></td><td><b>13</b></td></tr><tr><td>4.1 Elliptic-Hyperbolic approximation</td><td>13</td></tr></table>

<footnote>\*Aix Marseille Univ, CNRS, IUSTI, UMR 7343, Marseille, France, sergey.gavrilyuk@univ-amu.fr</footnote>

<footnote>†INRIA, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, 200 Avenue de la Vieille Tour, 33405 Talence cedex, France, mario.ricchiuto@inria.fr</footnote>
<table><tbody><tr><td>4.2</td><td>Hyperbolic relaxation</td><td>14</td></tr><tr><td>4.3</td><td>Numerical discretization</td><td>15</td></tr><tr><td>4.3.1</td><td>Time stepping</td><td>15</td></tr><tr><td>4.3.2</td><td>Discretization of the spatial differential operators</td><td>16</td></tr><tr><td><b>5</b></td><td><b>Numerical tests</b></td><td><b>17</b></td></tr><tr><td>5.1</td><td>Verification using analytical solutions</td><td>17</td></tr><tr><td>5.1.1</td><td>Solitary waves</td><td>17</td></tr><tr><td>5.1.2</td><td>Composite solution</td><td>17</td></tr><tr><td>5.2</td><td>Breakdown of a Gaussian water column</td><td>19</td></tr><tr><td>5.2.1</td><td>Validation against section averaged 2D shallow water simulations</td><td>20</td></tr><tr><td>5.3</td><td>Riemann problem</td><td>23</td></tr><tr><td>5.3.1</td><td>Validation against section averaged 2D shallow water simulations</td><td>24</td></tr><tr><td>5.3.2</td><td>Solitary wave fission</td><td>24</td></tr><tr><td>5.4</td><td>Treske's experiments in a trapezoidal channel</td><td>24</td></tr><tr><td><b>6</b></td><td><b>Conclusion</b></td><td><b>27</b></td></tr><tr><td><b>A</b></td><td><b>Geometrical coefficient properties for symmetric channels</b></td><td><b>28</b></td></tr><tr><td><b>B</b></td><td><b>Potential for symmetric channels</b></td><td><b>28</b></td></tr></tbody></table>

# 1 Introduction

This paper presents a new fully nonlinear one-dimensional dispersive model for gravity waves over topography in channels with prismatic sections. The model accounts for geometrical dispersion processes appearing in hydrostatic flows, and originated solely by bathymetric variations.

The physical background of this work is the propagation of bores in prismatic channels, and can be related to the laboratory experiments by [41], as well as to the in situ measurements in the Garonne and Seine rivers, reported in [7, 8]. Before the Froude number is high enough for wave breaking to occur (~1.3), two different dispersive propagation regimes are clearly observable both in laboratory and river experiments. One of them, occurring for Froude numbers higher than ~1.15, corresponds to the usual dispersive undular bore known since the experiments by [17]. These waves are well understood, and can be related theoretically to the classical dispersion relation from Airy theory using the analogy by [28]. For lower Froude numbers, longer waves are observed which do not correspond to classical dispersive propagation. This phenomenon has been studied in [13, 12] where it has been shown numerically and theoretically that these low Froude number undular bores are generated by a purely geometrical process related to refraction in the transverse direction. These waves are unaffected by vertical kinematics, and appear also in hydrostatic flows. So despite of their dispersive nature, these waves can be simulated using the hyperbolic shallow water equations as shown in [13, 12]. The dispersive behaviour of waves in materials with microstructure is well known, see e.g. [5] and [6]. Following the asymptotic analysis of [26] for anisotropic media, in [13, 12] the authors introduce the idea of a scale separation between the characteristic transverse length ($l$, running across the section), and the characteristic longitudinal length ($L$, in the main propagation direction along the channel axis). Using the small parameter $\epsilon = l/L$ the authors derive a dispersive linear wave equation providing an asymptotic section averaged approximation of the linearized shallow water system. The dispersion effects, and the dispersion coefficient, only depend on the channel's geometry. Using the theory by [28], the authors use the dispersion relation of their linear model to obtain a good quantitative prediction of the wavelengths of the low Froude number waves. These waves are purely related to geometrical effects, and in particular to bathymetric variations. They have been baptised "dispersive-like" in [12]. Later on, [33] used the same formal development to identify solitary waves associated to geometrical dispersion. They also introduce ad-hoc non-linear terms to obtain heuristic Boussinesq and KdV-like models.
From the modelling point of view, the first rigorous Boussinesq model for channels with non-rectangular sections dates to the early work by [32], focusing on prismatic channels. First generalizations allowing small section variations along the channel axis were proposed by [37, 38]. Equations accounting for arbitrary non-uniform cross sections have been developed in [43]. All the above models are obtained in a very general setting which starts from the Euler equations and then goes through classical scaling assumptions on length and velocity scales to arrive to models including weak dispersive effects. These scaling hypotheses differentiate between horizontal and vertical motions, and do not account for the differences between transverse and longitudinal horizontal processes. Despite of this general setting, it has been shown in [24] that e.g. the model by [43] only accounts for waves associated to geometrical variations, as the linear model by [12]. There is thus some ambiguity on which type of dispersion is modelled by existing equations, and which range of applications they can be used for.

Here we propose a 1D fully nonlinear model for dispersive-like waves originated from 2D bathymetric variations. In this respect, there present work is inherently different from previous ones on section averaged dispersive models. There is no ambiguity on the dispersive process modelled, which is fully relying on horizontal kinematics. The new geometrical Green-Naghdi model (gGN) model proposed is derived from the nonlinear hyperbolic shallow water equations. We call this model the Green-Naghdi type (and not, as conventionally, the Serre-Green-Naghdi type) to emphasize that it is based on the model taking into account the effects of 2D bathymetry. The physical assumptions made are somewhat similar to those invoked in [12], albeit retaining the full hydrodynamic non-linearity. The new section averaged gGN model is thus fully non-linear, it has a Lagrangian structure, and, as a consequence, an exact energy conservation law. We study the propagation properties of the model, and show the existence of several types of exact travelling wave solutions including solitary waves and periodic waves. We propose two reformulations to solve the PDE system numerically: the first based a elliptic-hyperbolic decomposition [19], the second using hyperbolic relaxation approach [18]. We discuss the numerical approximation of both formulations, and provide a thorough investigation or the model, which is validated both using comparisons with transverse averages of full 2D shallow water simulations, and against the experiments by [41] in the low Froude number range.

The structure of the paper is the following. In section 2 we introduce the asymptotic hypothesis and discuss the formal derivation of the model. Section 3 is devoted to the study of some of its properties: energy conservation law, dispersion relation, travelling wave solutions. In section 4 we propose some re-formulations of the system more suitable for its numerical approximation, and discuss the numerical methods implemented to this end. In section 5 we validate the model by comparing it with averaged 2D shallow water computations, as well as with experiments in prismatic channels with a trapezoidal cross-section. Technical details are given in Appendices.

## 2 Equations averaged in $y$-direction

We consider shallow water flows over topography described by the classical 2D Saint-Venant equations :

$$h_t + (hu)_x + (hv)_y = 0, \qquad (1a)$$

$$(hu)_t + (hu^2 + gh^2/2)_x + (huv)_y = -ghb_x, \qquad (1b)$$

$$(hv)_t + (huv)_x + (hu^2 + gh^2/2)_y = -ghb_y. \qquad (1c)$$

Here $t$ is time, $(x, y)$ are the Cartesian coordinates respectively running in the longitudinal (channel axis, $x$), and transverse (across the sections, $y$) directions (see the sketch on figure 1). We denote by $(u, v)$ the $x$ and $y$ velocity components, $b(x, y)$ is the bottom topography, and $g$ is the acceleration of gravity.

We assume that the transverse and longitudinal propagation processes are characterized by two different length scales $l$ and $L$, respectively, such that $l \ll L$. We thus introduce the small paramter

$$\varepsilon := \frac{l}{L}. \qquad (2)$$

Using the scaling (dimensionless variables are denoted with the ‘tilde’ sign)

$$x = L\tilde{x}, \quad y = \varepsilon L \tilde{y}, \quad t = L/\sqrt{gH} \tilde{t}, \quad u = \sqrt{gH} \tilde{u}, \quad v = \varepsilon \sqrt{gH} \tilde{v}, \\ h = H \tilde{h}, \quad b = H \tilde{b},$$

where $H$ is the characteristic water depth, and suppressing for simplicity the ‘tilde’, we obtain the following dimensionless equations

$$h_t + (hu)_x + (hv)_y = 0, \qquad (3a)$$

$$u_t + uu_x + vu_y + (h+b)_x = 0, \qquad (3b)$$

$$\varepsilon^2(v_t + uv_x + vv_y) + (h+b)_y = 0, \qquad (3c)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: Top: geometrical sketch of a channel with longitudinal (x), transverse (y) and vertical (z) axes. Bottom: typical channel sections of width l. Left: triangular section with $l = l_1 + l_2$, and $b(0) = b(l) = b_0$, $b(l_1) = 0$. Right: trapezoidal section with $l = l_1 + l_2 + l_3$ and $b(0) = b(l) = b_0$, $b(l_1) = b(l_1 + l_2) = 0$.</figcaption></figure>

We now introduce some of our main hypotheses. The first is that the channel has a prismatic bathymetry defined as $b = b(y)$. A slow variation of b with x could also be included, with for example $\partial_x b = O(\varepsilon^\beta)$ for some $\beta > 0$. We will leave this out for simplicity.

The second important hypothesis is that we assume $(hv)|_{y=0} = (hv)|_{y=l}$. This hypothesis is relevant for channels with banks ($h=0$ on both sides), as well as channels with reflective straight walls ($v=0$ on both sides), and of course the periodic case with also an $l$- periodic bathymetry $b(y) = b(l+y)$. Our aim is to propose equivalent 1D equations, which we will show exhibit a dispersive behaviour. To this end we will average (3) in the transverse direction. We will average across a length $l$ which we assume to be constant. This length can be a given channel width (for channels with vertical walls), or the width of a reference/initial wet section. There is no approximation for channels with walls, or periodic channels. Figure 1 shows typical geometries encountered in practice, and gives some initial notation.

## 2.1 Generalities and asymptotic expansion

For a given function $f(t, x, y)$ we define two types of averaging operators in $y$-direction. One is standard :

$$ \bar{f}(t, x) = \frac{1}{l} \int_{0}^{l} f(t, x, y) dy \quad (4) $$

The second one defines the weighted averaged value

$$ \langle f \rangle = \frac{1}{l\bar{h}} \int_{0}^{l} h(t, x, y) f(t, x, y) dy. \quad (5) $$

Sometimes the last averaging is called *Favre averaging*. The advantage of the last definition is that the mass conservation law is exactly satisfied in an usual sense, with the velocity $\langle u \rangle$. The dimensionless averaged conservative form of equations (1a), (1b) in the case $b = b(y)$ becomes :

$$ \bar{h}_t + (\bar{h}\langle u \rangle)_x = 0, \quad (6a) $$

$$ (\bar{h}\langle u \rangle)_t + (\overline{hu^2} + \frac{1}{2}\overline{h^2})_x = 0, \quad (6b) $$

We need to find $\overline{hu^2}$ and $\overline{h^2}$. Equation (3c) implies :

$$ h + b = \bar{h} + \bar{b} + O(\varepsilon^2). \quad (7) $$

We want to find this $O(\varepsilon^2)$ term in (7) in explicit form. For this, we make the following ansatz for the velocity field :

$$ u(t, x, y) = \langle u \rangle(t, x) + O(\varepsilon^2), \quad v = v_0(t, x, y) + O(\varepsilon). \quad (8) $$
This ansatz means that the shear effects in y-direction are negligible, i.e. the velocity in x-direction is almost y - independent (see [30] for the derivation of a model for fluid flows in prismatic channels with shear effects, but without dispersive effects). The anzatz (8) can be considered as a potential flow assumption commonly used for the derivation of the Saint-Venant equations. Indeed, the equations (3a)–(3c) imply the dimensionless vorticity equation in the form

$$ \frac{D}{Dt} \left( \frac{\omega}{h} \right) = 0, \quad \frac{D}{Dt} = \frac{\partial}{\partial t} + u \frac{\partial}{\partial x} + v \frac{\partial}{\partial y}, \quad \omega = \varepsilon^2 v_x - u_y. $$

If $\omega = 0$ initially, it stays zero for any time. Integrating the equation $\omega = 0$ with respect to $y$, we have

$$ u(t, x, y) = c(t, x) + \varepsilon^2 \int_0^y v_x(t, x, s) ds, $$

where $c(t, x)$ is a function of $(t, x)$. Multiplying the expression for $u$ by $h$, one obtains

$$ h(t, x, y) u(t, x, y) = h(t, x, y) c(t, x) + \varepsilon^2 h(t, x, y) \int_0^y v_x(t, x, s) ds. $$

The averaging with respect to $y$ gives then

$$ c(t, x) = \langle u \rangle(t, x) + \mathcal{O}(\varepsilon^2), $$

and

$$ u - \langle u \rangle = \mathcal{O}(\varepsilon^2). $$

Also,

$$ \begin{aligned} \overline{hu^2} &= \overline{h(u - \langle u \rangle + \langle u \rangle)^2} = \overline{h(u - \langle u \rangle)^2} + 2\overline{h(u - \langle u \rangle)} \langle u \rangle + \overline{h\langle u \rangle^2} \\ &= \overline{h(u - \langle u \rangle)^2} + \overline{h\langle u \rangle^2} = \overline{h\langle u \rangle^2} + \mathcal{O}(\varepsilon^4). \end{aligned} \quad (9) $$

An additional hypothesis about the behavior of $v_0$ will be now specified. Replacing ansatz (8) into the mass conservation law (3a), one obtains for $v_0$ the equation :

$$ \overline{h}_t + (\langle u \rangle (\overline{h} + \overline{b} - b))_x + (v_0 (\overline{h} + \overline{b} - b))_y = \mathcal{O}(\varepsilon). \quad (10) $$

Since

$$ \overline{h}_t + (\overline{h}\langle u \rangle)_x = 0, $$

one gets

$$ (\langle u \rangle (\overline{b} - b))_x + (v_0 (\overline{h} + \overline{b} - b))_y = \mathcal{O}(\varepsilon). $$

Neglecting the $\mathcal{O}(\varepsilon)$ term, one obtains the following zero order approximation :

$$ v_0(\overline{h} + \overline{b} - b) \approx C(t, x) - S(y) \langle u \rangle_x, \quad S(y) = \int_0^y (\overline{b} - b(y')) dy'. $$

The function $S(y)$ plays an important role in this work. It has the property $S(0) = S(l) = 0$. Moreover we prove in appendix A that in the case $b(y) = b(l - y)$ for any $0 < y < l$, i.e. $b(y)$ symmetric with respect to the channel axis, we have $\overline{S} = 0$. We will see the importance of this result shortly.

Concerning the integration ‘constant’ $C(t, x)$, we need to determine it using additional hypotheses. For the applications we have in mind, we can use the requirement of a vanishing average discharge in the y- direction :

$$ \int_0^l v_0 h dy = \int_0^l v_0 (\overline{h} + \overline{b} - b) dy = 0, $$

which leads to $C(t, x) = \overline{S} \langle u \rangle_x$.

## 2.2 Small-scale geometrical effects on the section average flow

Using the above results we can now write

$$ v_0(\overline{h} + \overline{b} - b) \approx (\overline{S} - S(y)) \langle u \rangle_x. \quad (11) $$

In particular, (11) and mass conservation (6a) give for the transverse velocity

$$ v_0(t, x, y) \approx \frac{(\overline{S} - S(y))}{\overline{h} + S_y} \langle u \rangle_x = \frac{(\overline{S} - S(y))\dot{\tau}}{1 + \tau S_y}, \quad (12) $$
with $\tau = 1/\bar{h}$, and the “dot” denotes total derivative along the longitudinal direction :

$$ \dot{\tau} = \frac{\partial \tau}{\partial t} + \langle u \rangle \frac{\partial \tau}{\partial x}. \qquad (13) $$

We set

$$ \sigma(y) = \bar{S} - S(y). \qquad (14) $$

The function $\sigma(y)$ has the following properties :

$$ \sigma(0) = \sigma(l) = \bar{S}, \quad \bar{\sigma} = 0, \quad \frac{d\sigma}{dy} = b(y) - \bar{b}, \quad \frac{d\bar{\sigma}}{dy} = 0. \qquad (15) $$

We can thus write the following $O(\varepsilon)$ approximation for $v_0$ :

$$ v_0(t, x, y) \approx \frac{\sigma(y)\dot{\tau}(t, x)}{1 - \tau(t, x)\frac{d\sigma(y)}{dy}}. \qquad (16) $$

We are now able to find the $O(\varepsilon^2)$ terms from (3c):

$$ (h + b)_y \approx -\varepsilon^2(v_t + uv_x + vv_y) \approx -\varepsilon^2(\dot{v}_0 + v_0(v_0)_y). \qquad (17) $$

One has

$$ \dot{v}_0(t, x, y) = \frac{\sigma \ddot{\tau}}{1 - \tau \frac{d\sigma}{dy}} + \frac{\sigma \frac{d\sigma}{dy} \dot{\tau}^2}{\left(1 - \tau \frac{d\sigma}{dy}\right)^2}, \qquad (18) $$

and

$$ h + b \approx \bar{h} + \bar{b} - \varepsilon^2 (M - \overline{M}), \quad M(t, x, y) = \int_0^y \dot{v}_0(t, x, s) ds + \frac{v_0^2(t, x, y)}{2}. \qquad (19) $$

Or

$$ h \approx \bar{h} - \frac{d\sigma}{dy} - \varepsilon^2 (M - \overline{M}). \qquad (20) $$

Thus

$$ \overline{h^2} = \left(\overline{h} - \frac{d\sigma}{dy}\right)^2 - 2\varepsilon^2 \overline{\left(h - \frac{d\sigma}{dy}\right)} (M - \overline{M}) + O(\varepsilon^4) \qquad (21) $$

$$ = \overline{h^2} + 2\varepsilon^2 M \overline{\frac{d\sigma}{dy}} + \text{const} + O(\varepsilon^4), \qquad (22) $$

because $\overline{\frac{d\sigma}{dy}} = 0$. From the above we can deduce that within $O(\varepsilon^4)$

$$ \left(\overline{\frac{h^2}{2}}\right)_x = \left(\overline{\frac{h^2}{2}}\right)_x + \varepsilon^2 \left(\overline{M \frac{d\sigma}{dy}}\right)_x, \qquad (23) $$

where the general expression for $M$ is

$$ M = \int_0^y \left( \frac{\sigma(s) \ddot{\tau}}{1 - \tau \frac{d\sigma(s)}{ds}} + \frac{\sigma(s) \frac{d\sigma(s)}{ds} \dot{\tau}^2}{\left(1 - \tau \frac{d\sigma}{ds}\right)^2} \right) ds + \frac{\sigma^2(y) \dot{\tau}^2}{2 \left(1 - \tau \frac{d\sigma}{dy}\right)^2}. \qquad (24) $$

Note that the denominator $1 - \tau \frac{d\sigma}{dy} = 1 - (b(y) - \bar{b})/\bar{h}$ can be easily shown to be strictly positive for many sections of practical relevance as e.g. trapezoidal ones.

We introduce two important particular cases, and then will devote the paper in-depth study of the one of them. We first consider channels verifying the assumption $\bar{S} = 0$ which is, in particular, satisfied for symmetric channel cross-sections (see appendix A for a proof). Formally, since only $O(\varepsilon^2)$ terms are retained in the equations, this assumption can be weakened to $\bar{S} = O(\varepsilon^\beta)$, with $\beta > 0$. In this case one could speak of quasi-symmetric channel cross-sections. This first family of models can be characterized by the following property.
**Proposition 1** (Lagrangian structure for symmetric and quasi-symmetric channel cross-sections). Consider quasi-symmetric channel cross-sections for which $\bar{S} = O(\varepsilon^\beta)$ with $\beta > 0$. Then, up to terms of order $O(\varepsilon^\beta)$, the term $\overline{M \frac{d\sigma}{dy}}$ can be written as the variational derivative of the Lagrangian potential

$$ \mathcal{L}(\tau, \dot{\tau}) = \overline{N \frac{d\sigma}{dy}} $$

with

$$ N(\tau, \dot{\tau}, y) = \frac{\dot{\tau}^2}{2} \int_0^y \frac{\sigma(s)ds}{1 - \tau \frac{d\sigma}{ds}}. \qquad (25) $$

More precisely, we have

$$ \overline{M \frac{d\sigma}{dy}} = - \left( \frac{\partial \mathcal{L}}{\partial \tau} - \frac{D}{Dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{\tau}} \right) \right) + \mathcal{O}(\varepsilon^\beta) \qquad (26) $$

For symmetric channel cross-sections verifying the condition $\bar{S} = 0$, the expression (26) is exact (the term $O(\varepsilon^\beta)$ identically vanishes).

The proof is given in appendix B.

The last proposition shows that for quasi-symmetric channel cross-sections the pressure term in our model has a Lagrangian structure :

$$ \frac{\overline{h^2}}{2} \approx \frac{(\overline{h})^2}{2} - \varepsilon^2 \frac{\delta \mathcal{L}}{\delta \tau}, \quad \text{with} \quad \frac{\delta \mathcal{L}}{\delta \tau} = \frac{\partial \mathcal{L}}{\partial \tau} - \frac{D}{Dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{\tau}} \right). \qquad (27) $$

We call this model symmetric geometrical Green-Naghdi equations. The advantage of the formulation (27) are multiple. In particular, the form (27) guarantees the variational formulation of our equations ([20]), and, as a consequence, the existence of an exact energy balance law. More details are given in the next Section.

The derivation above provides a large family of models. Its form however does not allow to obtain closed form expressions for the dispersion coefficients. To obtain some simplification we consider a different geometrical ansatz.

As we have already mentioned above, the condition $1 - \tau \frac{d\sigma}{dy} > 0$ can be easily shown for sections of practical importance. An interesting limit is obtained when considering channels which are wide compared to the depth. This typically the case in realistic applications both in natural and artificial environnements. For this case, we can assume further that ratio $\tau \frac{d\sigma}{dy}$ is small, namely $\tau \frac{d\sigma}{dy} = O(\varepsilon^\gamma)$ for some $\gamma > 0$. In this case one can easily show that

$$ M = \int_{0}^{y} \left( \sigma(s) \ddot{\tau} + \sigma(s) \frac{d\sigma(s)}{ds} \dot{\tau}^{2} \right) ds + \frac{\sigma^{2}(y)}{2} \dot{\tau}^{2} + \mathcal{O}(\varepsilon^{\gamma}) $$

Neglecting the asymptotically small $O(\varepsilon^\gamma)$ term, we obtain

$$ M \approx \ddot{\tau} \int_{0}^{y} \sigma(s)ds + \sigma^{2}(y)\dot{\tau}^{2} - \dot{\tau}^{2} \frac{\sigma^{2}(0)}{2}. $$

Hence,

$$ \overline{\frac{d\sigma}{dy} M} = \ddot{\tau} \overline{\frac{d\sigma}{dy} \int_{0}^{y} \sigma(s) ds} = \ddot{\tau} \left( \overline{\frac{d}{dy} \left( \sigma(y) \int_{0}^{y} \sigma(s) ds \right)} - \sigma^2 \right) \qquad (28) $$

$$ = \ddot{\tau} \left( \sigma(l)\bar{\sigma} - \overline{\sigma^2} \right) = -\ddot{\tau}\sigma^2 = -\chi\ddot{\tau}, \quad \text{with} \quad \chi = \overline{S^2} - \overline{S^2} > 0. \qquad (29) $$

This limit also enjoys a very important structural property which is resumed in the following result.

**Proposition 2** (Lagrangian structure for wide channel cross-sections). Consider wide channels for which $\tau \frac{d\sigma}{dy} = O(\varepsilon^\gamma)$, with $\gamma > 0$. Then, up to terms of order $O(\varepsilon^\gamma)$, the hydrodynamic term $\overline{M \frac{d\sigma}{dy}}$ can be written as the variational derivative of a Lagrangian potential. More precisely we have

$$ \overline{M \frac{d\sigma}{dy}} = - \left( \frac{\partial \mathcal{L}}{\partial \tau} - \frac{D}{Dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{\tau}} \right) \right) + \mathcal{O}(\varepsilon^\gamma) \qquad (30) $$
with

$$ \mathcal{L}(\tau, \dot{\tau}) = -\chi \frac{\dot{\tau}^2}{2}. \qquad (31) $$

and $\chi$ as in (29).

*Proof.* The proof follows is directly from the explicit expression (31) of $\mathcal{L}(\tau, \dot{\tau})$. $\square$

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

We will refer to the last model, having a very simple form. as simplified geometrical Green-Naghdi equations. We will devote the remainder of the paper to the study of its theoretical properties, and to its applications to the propagation of undular bores in channels.

The simplified geometrical Green-Naghdi model can be considered as a prototype averaged model, containing only one (geometrical) parameter $\chi$. In the case of symmetric channel cross-sections, $\chi = \overline{S^2}$.

For the prototype average model equations (6a) - (6b) can be written in a simple closed form reading, in dimensional variables :

$$ \bar{h}_t + (\bar{h}\langle u \rangle)_x = 0, \qquad (32) $$

$$ (\bar{h}\langle u \rangle)_t + (\bar{h}\langle u \rangle^2 + g\frac{\bar{h}^2}{2} + p)_x = 0, \qquad (33) $$

where we set

$$ p = -\chi\ddot{\tau}, \quad \tau = \frac{1}{\bar{h}}, \quad \chi = \overline{S^2} - \overline{S}^2 > 0. \qquad (34) $$

with $\chi$ having dimensions [$m^4$].

## 2.3 Expressions for the geometric parameter for the triangular and trapezoidal sections

The value of the parameter $\chi$ for a triangular type bottom (see Figure 1) is :

$$ \chi = \frac{b_0^2(l_1^2 + 4l_1l_2 + l_2^2)}{720}. \qquad (35) $$

Here $b_0$ is the bottom height ($b_0 = b(0) = b(l)$), and $l_1$ is the position of the bottom singularity ($b(l_1) = 0$), and $l = l_1 + l_2$). The value of the parameter $\chi$ for a trapezoidal type bottom is :

$$ \chi = \frac{b_0^2(l_1 + l_3)(l_1^3 + l_1^2(6l_2 + 5l_3) + l_3(15l_2^2 + 6l_2l_3 + l_3^2) + l_1(15l_2^2 + 24l_2l_3 + 5l_3^2))}{720(l_1 + l_2 + l_3)^2}. \qquad (36) $$

Here $b_0$ is the bottom height ($b_0 = b(0) = b(l)$), $l_1$ is the position of the first bottom singularity ($b(l_1) = 0$), $l_1 + l_2$ is the position of the second bottom singularity ($b(l_1 + l_2) = 0$), $b(y) = 0$ for $y \in (l_1, l_1 + l_2)$, and $l = l_1 + l_2 + l_3$.

The formulas coincide if one takes $l_2 = 0$ and replaces $l_3$ by $l_2$. Moreover the two formulas are symmetric with respect the pairs $(l_1, l_2)$ and $(l_1, l_3)$, respectively. This implies that the dispersion effects are insensitive to mirroring of asymmetric shapes, as one should expect. Similarly, one easily verifies that asymmetric configurations have a smaller dispersion coefficient than symmetric ones. For this, one has to look for the minimum of $\chi$ under the constraint $l_1 + l_2 = l$ in the case of the triangular type bottom, or under the constraint $l_1 + l_2 + l_3 = l$ in the case of trapezoidal type bottom. For the triangular bottom topography, the maximal value of $\chi/(b_0^2l^2)$ is achieved for $l_1 = l_2 = l/2$, and is equal approximately to 0.0020833. For the trapezoidal bottom the maximal value of $\chi/b_0^2l^2$ is achieved for $l_1 = l_3 = (21 - \sqrt{41})l/40$, and is equal approximately to 0.0027139.

# 3 Properties of the geometrical Green-Naghdi type system

## 3.1 Galilean invariance

The geometrical models obtained in the previous section can be formally written as

$$ \dot{\bar{h}} + \bar{h}\langle u \rangle_x = 0, $$

$$ \bar{h}\langle \dot{u} \rangle + g\bar{h}\bar{h}_x + p(\tau, \dot{\tau}, \ddot{\tau})_x = 0, $$

These equations involve only space derivative and first and second 1D material time derivatives (denoted, as usually, by ‘dot’). Now, consider a Galilean reference frame moving with a constant velocity. Both space and material time derivatives are invariant wrt this change of the reference frame. Thus the model is Galilean invariant.
## 3.2 Lagrangian structure and energy equation

Consider the Hamilton action

$$a = \int_{t_0}^{t_1} L \, dt, \quad L = \int_{-\infty}^{+\infty} \left( \bar{h} \frac{\langle u \rangle^2}{2} - W(\bar{h}, \dot{\bar{h}}) \right) dx, \quad W(\bar{h}, \dot{\bar{h}}) = \bar{h} w(\bar{h}, \dot{\bar{h}}). \quad (37)$$

Here $t_i, i = 0, 1$, are given time instants, and $w(\bar{h}, \dot{\bar{h}})$ is a given potential. The mass conservation law

$$\bar{h}_t + (\bar{h}\langle u \rangle)_x = 0$$

is considered as a differential constraint. Then, the Euler-Lagrange equation (momentum equation) for (37) can be written as (see e.g. [20] and [39]):

$$(\bar{h}\langle u \rangle)_t + \left(\bar{h}\langle u \rangle^2 + \bar{h}^2 \frac{\delta w}{\delta \bar{h}}\right)_x = 0, \quad \text{with} \quad \frac{\delta w}{\delta \bar{h}} = \frac{\partial w}{\partial \bar{h}} - \frac{D}{Dt} \left(\frac{\partial w}{\partial \dot{\bar{h}}}\right).$$

Introducing a new potential $\tilde{w}(\tau, \dot{\tau})$ as:

$$\tilde{w}(\tau, \dot{\tau}) = \tilde{w}\left(\frac{1}{\bar{h}}, -\frac{\dot{\bar{h}}}{\bar{h}^2}\right) = w(\bar{h}, \dot{\bar{h}}),$$

one can also rewrite the momentum equation in the form

$$(\bar{h}\langle u \rangle)_t + \left(\bar{h}\langle u \rangle^2 - \frac{\delta \tilde{w}}{\delta \tau}\right)_x = 0, \quad \text{with} \quad \frac{\delta \tilde{w}}{\delta \tau} = \frac{\partial \tilde{w}}{\partial \tau} - \frac{D}{Dt} \left(\frac{\partial \tilde{w}}{\partial \dot{\tau}}\right).$$

Thus, the averaged equations in the transverse direction both for the symmetrical case and the simplified case are obtained by taking

$$\tilde{w}(\tau, \dot{\tau}) = \frac{g}{2\tau} + \mathcal{L}(\tau, \dot{\tau}).$$

Setting $\bar{q} = \bar{h}\langle u \rangle$, the system admits the energy conservation law which is a consequence of invariance of the Lagrangian under time shift:

$$\left(\bar{h}\left(\frac{\langle u \rangle^2}{2} + g\frac{\bar{h}}{2} + E\right)\right)_t + \left(\bar{q}\left(\frac{\langle u \rangle^2}{2} + g\frac{\bar{h}}{2} + E\right) + \left(g\frac{\bar{h}^2}{2} - \frac{\delta \mathcal{L}}{\delta \tau}\right)\langle u \rangle\right)_x = 0, \quad (38)$$

with

$$E = \mathcal{L} - \dot{\tau} \frac{\partial \mathcal{L}}{\partial \dot{\tau}}.$$

In particular, for model (32)-(33)-(34) we have

$$\left(\bar{h}\left(\frac{\langle u \rangle^2}{2} + g\frac{\bar{h}}{2} + \frac{\chi \dot{\tau}^2}{2}\right)\right)_t + \left(\bar{q}\left(\frac{\langle u \rangle^2}{2} + g\frac{\bar{h}}{2} + \frac{\chi \dot{\tau}^2}{2}\right) + \left(g\frac{\bar{h}^2}{2} - \chi \ddot{\tau}\right)\langle u \rangle\right)_x = 0. \quad (39)$$

## 3.3 Linear dispersion properties

Consider now the linearized version of (32), (33) and (34):

$$\bar{h}_t + \bar{h}_0 \langle u \rangle_x = 0, \quad (40)$$

$$\langle u \rangle_t + g \bar{h}_x = \frac{\chi}{\bar{h}_0^2} \langle u \rangle_{xxt}. \quad (41)$$

Using standard techniques we show that the dispersion relation of the system is

$$\omega^2 = k^2 \frac{c_0^2}{1 + \frac{\chi}{\bar{h}_0^2} k^2} \quad (42)$$

with $k$ the wavenumber and $c_0 = \sqrt{g\bar{h}_0}$. In figure 2 we show the typical behaviour of the non-dimensional phase speed $\omega/kc_0$ as a function of the horizontal reduced wave number $kl$, for the case of a trapezoidal section. We
<figure><img src="image_4.png" /><figcaption>Figure 2: Dispersion relation for trapezoidal sections. Left: notation. Center: deep channel with steep slopes. Right: shallow channel with mild slopes. Comparison of the simplified gGN model (in blue) with the phase celerity of the linear model of [12] (in red).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Lemoine analogy for the experiments by [41]. Comparison of the theoretical predictions using the dispersion relation of the simplified gGN model (blue), with the one of [12] (red), and the data by [41] (symbols).</figcaption></figure>

consider the case of a deep channel with steep slopes, and of a shallow channel with mild slope. Using the notation of the leftmost picture in the figure, these cases correspond respectively to $(\beta, w, b_0) = (1/2, 2m, 2m)$, and $(\beta, w, b_0) = (1/3, 1.24m, 0.2m)$. The latter is representative of the channel used in the well known experiments by [41]. The coefficient $\chi$ is computed from (36) with $l_2 = w$, $l_1 = l_3 = b_0/\beta$. For completeness, we compare the curves with those obtained using the dispersion relation of the linear model obtained in [12]. The curves show comparable trends, with some differences for shorter waves which are within classical modelling error. In particular, for reduced wave numbers below $kl = \pi$ are within 2% in the first case, and 0.5% for the mild slopes.

For the second case, representative of the experiments by [41], we can use the dispersion relation to predict the wave length-Froude number dependence. This is achieved by equating for each Froude number the phase speed to the relative bore. We refer to [12, 28] for details. The wave lengths (scaled by the pre-bore section averaged depth $\bar{h}_1$) predicted theoretically using this method are reported on figure 3. For completeness we also report in the figure the data by [41], which include both those associated to geometrical dispersion (longer waves), and those measured on the channel axis for high Froude, related to classical non-hydrostatic dispersive effects. The result of figure 3 shows that the new model provides a good prediction of the long waves associated to geometrical dispersion. The differences wrt the linear model of [12] for higher Froude numbers is in line with the fact that the simplified gGN provides slightly larger phase celerities for shorter waves, as seen in figure 2. Nonlinear simulations of these experiments are discussed in section §5.4.

### 3.4 Travelling wave solutions: solitary waves

We show in this section that the system derived admits exact travelling wave solutions. We focus here on solitary waves. We consider solitary wave of amplitude $a$ and celerity $c$ travelling on a far field with constant averaged depth $\bar{h}_\infty$. So, let us set $\xi = x - ct$ and assume that at time $t = 0$ the peak is at $x = 0$. We now seek solutions of the form

$$ \bar{h} = \bar{h}(\xi), \quad \langle u \rangle = \langle u \rangle(\xi) $$
We can readily write from the mass conservation law that

$$ \bar{h}(\langle u \rangle - c) = m = \text{const}, $$

Which using the conditions at infinity gives the relations

$$ \langle u \rangle = \left(1 - \frac{\bar{h}_{\infty}}{\bar{h}}\right)c, \quad m = -\bar{h}_{\infty}c \qquad (43) $$

The momentum equation implies :

$$ \bar{h}(\langle u \rangle - c)^2 + \frac{g\bar{h}^2}{2} - \chi\ddot{\tau} = q = \text{const.} \qquad (44) $$

Since $\ddot{\tau} = (\langle u \rangle - c)((\langle u \rangle - c)\tau')' = m^2\tau(\tau\tau)'$ (here 'prime' means the derivative with respect to $\xi$), one obtains :

$$ m^2\tau + \frac{g}{2\tau^2} - \chi m^2\tau(\tau\tau)' = q = \text{const.} \qquad (45) $$

Let $\tau_{\infty} = 1/\bar{h}_{\infty}$. Then one can find $q$ and we have :

$$ m^2(\tau - \tau_{\infty}) + \frac{g}{2} \left( \frac{1}{\tau^2} - \frac{1}{\tau_{\infty}^2} \right) - \chi m^2 \tau (\tau \tau)' = 0. \qquad (46) $$

Multiplying by $\tau'$, one can integrate it again :

$$ m^2(\tau - \tau_{\infty})^2 + g \left( -\frac{1}{\tau} + \frac{1}{\tau_{\infty}} - \frac{\tau}{\tau_{\infty}^2} + \frac{1}{\tau_{\infty}} \right) - \chi m^2 (\tau \tau')^2 = 0. $$

We rewrite this equation as

$$ \chi m^2 \tau^2 \tau'^2 = F(\tau), \quad F(\tau) = m^2(\tau - \tau_{\infty})^2 - \frac{g}{\tau \tau_{\infty}^2} (\tau - \tau_{\infty})^2. \qquad (47) $$

On can remark that that at $\tau = \tau_{\infty}$, one has $F = 0$, $F' = 0$. Moreover, $F'' > 0$ if $m^2 > g/\tau_{\infty}^3$. If the velocity $\langle u \rangle_{\infty} = 0$, the last inequality means that the Froude number should be large than one (a classical inequality for the existence of solitary waves). Since $F \to -\infty$ as $\tau \to 0$, it implies that there exists $\tau_0 < \tau_{\infty}$ such that $F(\tau_0) = 0$. It implies the existence of supercritical solitary waves having the amplitude $a = \bar{h}_0 - \bar{h}_{\infty}$, $\bar{h}_0 = 1/\tau_0$.

We can now evaluate the above relation at the peak, where $\tau' = 0$ we have $\tau_0 = 1/\bar{h}_0 = 1/(\bar{h}_{\infty} + a)$. We thus obtain a relation celerity/amplitude:

$$ 0 = m^2 \frac{a^2}{\bar{h}_{\infty}^2 (\bar{h}_{\infty} + a)^2} - \frac{ga^2}{\bar{h}_{\infty} + a} $$

Using the constancy of $m$ we can also deduce $m^2 = c^2\bar{h}_{\infty}^2$ which allows to derive the compatibility relations

$$ c^2 = g(\bar{h}_{\infty} + a) = g\bar{h}_0 = \frac{g}{\tau_0}, \quad m^2 = g(\bar{h}_{\infty} + a)\bar{h}_{\infty}^2 = \frac{g}{\tau_0\tau_{\infty}^2} \qquad (48) $$

Combining all of the above results we end with the nonlinear ODE

$$ \tau'^2 = \frac{(\tau - \tau_{\infty})^2 (\tau - \tau_0)}{\chi \tau^3}. \qquad (49) $$

The simplest way to solve it is to construct only a half-width solitary wave with initial condition $\tau = \tau_0$ at $\xi = 0$.

In practice, we have actually solved (46) using the basic odeint method in numpy, starting from the peak with given zero derivative, and the known value of $\tau = \tau_0$. Typical shapes of the solitary waves obtained for two different values of the non-linearity $\epsilon = a\tau_{\infty}$, and of the dispersion parameter $\chi$ are reported in figure 4. For sake of comparison we report the SGN solitary wave with the same amplitude. The solitary wave of the new system generally have a more “peaky” shape, and a strong dependence on the relative values of $\chi$ and $\epsilon$.
<figure><img src="image_6.png" /><figcaption>Figure 4: Solitary waves obtained for different values of $\chi$, and for $\epsilon = 0.2$ (left) and $\epsilon = 0.5$ (right).</figcaption></figure>

## 3.5 Travelling waves solutions: a composite solution

Equation (45) allows to obtain solutions of other forms than solitary waves, including e.g. periodic solutions. Following the approach by [21], we exploit this property to show the existence of exact composite travelling waves involving two moving profiles connected via a generalized Rankine-Hugoniot relation. One of the two solutions is a periodic oscillation verifying a generalized form of equation (45) which can be written as

$$ \chi \tau^3(\tau')^2 = (\tau - \tau_1)(\tau - \tau_2)(\tau - \tau_3), \quad 0 < \tau_1 < \tau_2 < \tau_3, \quad m^2 \tau_1 \tau_2 \tau_3 = g. \qquad (50) $$

If $\tau_2 = \tau_3 = \tau_\infty$, one obtains (49). Otherwise equation (50) describes solutions oscillating between $\tau_1$ and $\tau_2$. To construct a stable solution we proceed as follows (see [21] for details)

1. Choose three roots such that $0 < \tau_1 < \tau_2 < \tau_3$, with $\tau_2$ and $\tau_3$ close to each other.

2. Choose a sign for the mass flow obtained from the second in (50). For $m < 0$ we obtain a right travelling wave.

3. Give an arbitrary value to the velocity corresponding to the state $\tau_2$. For example, $u_2 = 0$. Compute the celerity $C$ of the periodic traveling wave from $\bar{h}_2(u_2 - C) = m, \bar{h}_2 = 1/\tau_2$. Thus, the periodic traveling wave is perfectly defined.

4. From the first in (50) we obtain the value of the second derivative $\tau''|_{\tau=\tau_1}$:

$$ 2\chi\tau_1^3\tau''|_{\tau=\tau_1} = (\tau_1 - \tau_2)(\tau_1 - \tau_3) > 0. \qquad (51) $$

5. Using the generalized Rankine-Hugoniot relations (including dispersive terms), compute the state with index $\star$ (assumed constant in space):

$$ m = (u_{\star} - C)/\tau_{\star}, \qquad (52) $$

$$ m^2\tau_1 + \frac{g}{2\tau_1^2} - \chi m^2\tau_1^2\tau''|_{\tau=\tau_1} = m^2\tau_{\star} + \frac{g}{2\tau_{\star}^2}. \qquad (53) $$

Take the root $\tau_\star$ closest to $\tau_2$.

6. Then the solution having the constant state $\star$ on the left and the corresponding periodic traveling wave on the right, is stable: it contains a jump at which the generalized Rankine-Hugoniot conditions are satisfied.

7. Introduce a regularization of the jump using e.g. a half solitary wave from $h_\star = 1/\tau_\star$ to $h_1 = 1/\tau_1$ (in h variables).

To build the above solution in practice, given as data the tuple $(\tau_1, \tau_2, \tau_3)$, we first solve for $\tau_\star$, and then integrate (50) with initial value $\tau_1$ to compute the periodic branch. The integration is done in practice recasting (50) as

$$ 2\chi(\tau\tau')' = (1 - \frac{\tau_2}{\tau})(1 - \frac{\tau_3}{\tau}) + (1 - \frac{\tau_1}{\tau})(1 - \frac{\tau_3}{\tau}) \\ + (1 - \frac{\tau_1}{\tau})(1 - \frac{\tau_2}{\tau}) - (1 - \frac{\tau_1}{\tau})(1 - \frac{\tau_2}{\tau})(1 - \frac{\tau_3}{\tau}) \qquad (54) $$

which is integrated with initial conditions $(\tau, \tau') = (\tau_1, 0)$. The results presented have been obtained using the 8th order Runge-Kutta method in the SciPy class `solve_ivp`, with the relative and absolute tolerances `rtol=10⁻⁸` and `atol=10⁻¹²`. Concerning the regularization, we have computed the half solitary wave as discussed in the previous section.
<figure><img src="image_7.png" /><figcaption>Figure 5: Composite solution obtained for $\tau_1 = 1\text{m}^{-1}$, $\tau_2 = 1.3\text{m}^{-1}$, and $\tau_3 = 1.301\text{m}^{-1}$, and $\chi = 0.4\text{m}^4$. Left: entire solution computed on the domain [-350, 2100]m. Right: zoom of the regularized jump and first peak</figcaption></figure>

We show an example obtained choosing $(\tau_1, \tau_2, \tau_3) = (1\text{m}^{-1}, 1.3\text{m}^{-1}, 1.301\text{m}^{-1}$, and setting $\chi = 0.4\text{m}^4$. For these choices we obtain $\tau_* = 1.30050031985595\text{m}^{-1} < \tau_1$. Integrating the ODEs we obtain the solution plotted on figure 5. In the figure we also compare three treatments of the transition: the one obtained using a gGN solitary wave, a Gaussian $e^{-(x/L)^2}$ with $L = 8\text{m}$, and no regularization. We will compare the evolution of the three later in the results section.

# 4 System reformulations and numerical approximation

We propose in this section two re-formulations of the model in view of its numerical discretization. The first is based on a splitting of the original system in its elliptic and hyperbolic parts. This is a rather classical way to solve the equations, used in many works over the years [42, 27, 9, 31, 23]. More specifically, here we follow the reformulation proposed in [19, 25] for the Serre-Green-Naghdi equations, and in [11, 40, 24] for other types of Boussinesq models.

The second approach is based on a hyperbolization of the dispersive system which follows the early works by [29, 18, 10]. Since the homogenized equations (32), (33), (34) admit a variational formulation (37), the hyperbolic approximation can be derived by using the method of augmented Lagrangian proposed in [18, 14]. This method has been justified for the Serre-Green-Naghdi equations (the case of flat bottom) in [15]. Compared to classical relaxation methods, the method of augmented Lagrangian conserves the variational structure of the equations, thus guaranteeing the existence of an exact energy conservation law, as a direct consequence of the Noether theorem.

To simplify the notations, in this section we will simply denote $\bar{h}$ and $\langle u \rangle$ for the average depth and Favre averaged velocity, as $h$ and $u$.

## 4.1 Elliptic-Hyperbolic approximation

We consider a first rewriting of the system as a perturbation to the shallow water equations via an algebraic source. The compatibility with the original equations leads to an elliptic problem for the latter. To obtain this form, we start by using the relation

$$h^2 \ddot{\tau} = h(u_x)^2 + h(u_{xt} + uu_{xx}) = h(u_x)^2 + h\dot{u}_x = h(u_t + uu_x)_x = h(\dot{u})_x \quad (55)$$

where

$$\dot{u}_x = (u_{xt} + uu_{xx}) \neq (\dot{u})_x = (u_t + uu_x)_x.$$

Using the above, equation (33) can be recast as

$$(hu)_t + (hu^2 + g\frac{h^2}{2})_x = [\frac{\chi}{h}(\dot{u})_x]_x \quad (56)$$

Now following [19] we set

$$h\phi := (hu)_t + (hu^2 + gh^2/2)_x \Rightarrow \phi - ghx = \dot{u} \quad (57)$$

We can now combine the above definition with the compatibility requirement

$$h\phi = [\frac{\chi}{h}(\dot{u})_x]_x$$
and obtain the equation for $\phi$:

$$h\phi = \left[ \frac{\chi}{h} (\phi - gh_x)_x \right]_x.$$

The dispersive system (32), (33), (34) can thus be recast as

$$\begin{aligned} h\phi - \chi\left[\frac{1}{h}(\phi - gh_x)_x\right]_x &= 0 \\ h_t + (hu)_x &= 0 \\ (hu)_t + (hu^2 + gh^2/2)_x &= h\phi \end{aligned} \qquad (58)$$

which gives two operators. The first an elliptic problem for the auxiliary variable $\phi$, whose variational form reads

$$\int v h \phi + \chi \int v_x \tau \phi_x = \chi \int v_x \tau \delta_x \qquad (59)$$

where $v$ is a test function with compact support $v \in H_0^1$, and having set $\delta := gh_x$. We can have the following simple characterization.

**Proposition 3** (Coercivity of the elliptic operator). Provided there exist solution independent uniformly bounded values $h_{\min}$, and $h_{\max}$ such that $h \in (h_{\min}, h_{\max})$, then the variational form (59) is uniformly coercive.

*Proof.* Under the hypotheses made, we can easily show the $H^1$ norm equivalence

$$\frac{\chi C_{\min}}{h_{\max}} \|\phi\|_{H^1}^2 \le \int h\phi^2 + \chi \int \tau(\phi_x)^2 \le \frac{\chi C_{\max}}{h_{\min}} \|\phi\|_{H^1}^2$$

for any $C_{\max} \ge \max(1, \frac{h_{\min}h_{\max}}{\chi})$ and $C_{\max} \le \min(1, \frac{h_{\min}h_{\max}}{\chi})$. $\square$

The last two equations in (58) define an hyperbolic evolution operator for the shallow water variables, with an extra algebraic forcing term linear in $\phi$. This system can be effectively integrated in time and space.

## 4.2 Hyperbolic relaxation

In this method we introduce a second order equation for a new variable $\eta(t, x)$ which approximates the variable $\tau(t, x) = 1/h(t, x)$. The expression for $\dot{\eta}$ is written as a first order system if PDEs in the following way

$$\begin{aligned} (h\eta)_t + (hu\eta)_x &= h w, \\ (hw)_t + (huw)_x &= -\frac{\mu}{\chi}(h\eta - 1). \end{aligned}$$

Here $\mu$ is a large positive parameter, which allows to relax $\eta$ to $\tau(t, x)$. The full system now reads

$$\begin{aligned} h_t + (hu)_x &= 0, \\ (hu)_t + (hu^2 + p_{tot})_x &= 0, \\ (h\eta)_t + (hu\eta)_x &= hw, \\ (hw)_t + (huw)_x &= -\frac{\mu}{\chi}(h\eta - 1). \end{aligned} \qquad (60)$$

with the closure relation

$$p_{tot} = g \frac{h^2}{2} + \mu \tau (h\eta - 1) = g \frac{h^2}{2} + \mu (\eta - \tau). \qquad (61)$$

The system is hyperbolic with eigenvalues

$$\lambda_1 = u - c, \quad \lambda_{2,3} = u, \quad \lambda_4 = u + c, \qquad (62)$$

having set

$$c^2 = gh + \mu\tau^2. \qquad (63)$$

Following [14] and [39], one can easily show that (60) are the Euler-Lagrange equations for the Lagrangian :

$$L = \int_{-\infty}^{+\infty} \left( h \frac{u^2}{2} - W(h, \eta, \dot{\eta}) \right) dx,$$
with

$$W(h, \eta, \dot{\eta}) = \frac{gh^2}{2} - \frac{\chi h \dot{\eta}^2}{2} + \frac{\mu h}{2} \left( \eta - \frac{1}{h} \right)^2.$$

The pressure $p_{tot}$ is defined as the partial Legendre transform with respect to $h$:

$$p_{tot} = h \frac{\partial W}{\partial h} - W = \frac{gh^2}{2} + \mu (\eta - \tau), \qquad (64)$$

and the equation for $\eta$ coming from two last equations of (60):

$$\ddot{\eta} + \frac{\mu}{\chi}(\eta - \tau) = 0 \qquad (65)$$

is equivalent to the Euler-Lagrange equation for $\eta$:

$$\frac{\partial W}{\partial \eta} - \frac{\partial}{\partial t} \left( \frac{\partial W}{\partial \dot{\eta}} \right) - \frac{\partial}{\partial x} \left( \frac{\partial W}{\partial \dot{\eta}} u \right) = 0. \qquad (66)$$

In the limit $\mu \to \infty$ we recover the original system, and, in particular, the total pressure $p_{tot}$ becomes

$$p_{tot} = \frac{gh^2}{2} - \chi \ddot{\tau}. \qquad (67)$$

The interested reader can refer to [14, 39] and references therein for more details.

## 4.3 Numerical discretization

### 4.3.1 Time stepping

*Elliptic-Hyperbolic approximation.* For system (58) we start from initial conditions simply set on physical quantities $(h, u)$. The system is then evolved in time by means of an explicit multistage method. Setting

$$U_{EH} = (h, hu)^T, \qquad (68)$$

given the depth and velocities at time $n$ and at intermediate stages $(0, k)$, as well as $\phi^n, \dots, \phi^{k-1}$, the generic stage $k$ reads:

$$h^k \phi^k - \chi \left[ \frac{1}{h^k} (\phi^k - gh_x^k)_x \right]_x = 0 \quad \longrightarrow \quad \phi^k \\ U_{EH}^{k+1} = \sum_{\ell \ge 1} \alpha_\ell U_{EH}^{k+1-\ell} - \Delta t \sum_{\ell \ge 1} \beta_\ell \Psi_{SW}(U_{EH}^{k+1-\ell}, \phi^{k+1-\ell}), \qquad (69)$$

where $\Psi_{SW}$ is a discrete approximation of the shallow water spatial residual, including the extra source $\phi$. In this approach we need to assemble and invert the elliptic operator at each intermediate stage. In practice here we use (69) with two stages.

*Hyperbolic system.* In this case initial conditions are set on physical quantities $(h, u)$ as well as the auxiliary variables $w$ and $\eta$. For the latter initial conditions are deduced from the limit $\mu \to \infty$:

$$\eta(0, x) = \tau(0, x), \quad w(0, x) = (u(0, x))_x / h(0, x).$$

For solitary waves, relations (43) and (49) have been used to evaluate these quantities.

System (60) is then integrated in time using the method proposed in [18]. In particular, setting

$$U_{Hyp} = (h, hu, h\eta, hw)^T, \qquad (70)$$

given the solution value $U_{Hyp}^n$ at time $t^n$, we use a classical second order Strang splitting method:

$$\begin{align} U_{Hyp}^{n+1} &= O_1(\Delta t/2) O_{Hyp}(\Delta t) O_1(\Delta t/2) \nonumber \\ &\quad + O_2(\Delta t/4) O_{Hyp}(\Delta t/2) O_2(\Delta t/4) \nonumber \\ &\quad + \dots \nonumber \\ &\quad + O_5(\Delta t) O_{Hyp}(\Delta t) O_5(\Delta t), \tag{71} \end{align}$$

where the operator $O_1$ is the system of ODEs
Given initial data $(\tau_0, u_0, \eta_0, w_0)$ the above is integrated exactly as

$$ 
\begin{aligned}
\eta &= \eta_0 \cos(\sqrt{\alpha}t) + \tau_0 (1 - \cos(\sqrt{\alpha}t)) + \frac{w_0}{\sqrt{\alpha}} \sin(\sqrt{\alpha}t) \\
w &= w_0 \cos(\sqrt{\alpha}t) - (\eta_0 - \tau_0)\sqrt{\alpha} \sin(\sqrt{\alpha}t)
\end{aligned}
\qquad (72) $$

with $\alpha = \mu/\chi$.

The operator $O_{\text{Hyp}}$ is nothing else than the hyperbolic homogenous part of system (60), which is integrated using a multi-stage method (with the same notation of (69))

$$ U_{\text{Hyp}}^{k+1} = \sum_{\ell \ge 1} \alpha_{\ell} U_{\text{Hyp}}^{k+1-\ell} - \Delta t \sum_{\ell \ge 1} \beta_{\ell} \Psi_{\text{H}}(U_{\text{Hyp}}^{k+1-\ell}), \qquad (73) $$

where $\Psi_L$ is a discrete approximation of the right hand side of the homogenous hyperbolic system. In practice here we use a two stage method.

### 4.3.2 Discretization of the spatial differential operators

Spatial differential operators are approximated on a uniform discretization of the 1D spatial domain, with size $\Delta x = 1/N$, $N$ representing the number of mesh cells.

The elliptic operator in (69) is discretized with $P^1$ finite elements, using the variational formulation (59). The invertibility of the resulting problem is characterized by proposition 3. In practice, the system being tri-diagonal, we have use Thomas method which allows to obtain the solution in two explicit seeps. The third derivative of $h$ appearing in the right hand side is evaluated with a second order formula.

The hyperbolic operators, are discretized using a second order explicit Residual Distribution method [36, 1, 2]. To first order of accuracy in space one stage of the method can be written as

$$ U_i^{k+1} = U_i^n - \frac{\Delta t}{\Delta x} \phi_i^{i-1/2}(U^n) - \frac{\Delta t}{\Delta x} \phi_i^{i+1/2}(U^n) \qquad (74) $$

where the fluctuations $\Phi_i^{i\mp 1/2}(U^n)$ are defined by some upwind biased distribution of the full residuals in mesh cells $[x_{i-1}, x_i]$ and $[x_i, x_{i+1}]$. In particular, for a hyperbolic system reading

$$ U_t + F_x = S, \quad A = F_U $$

we define

$$ \phi^{i-1/2} = \int_{x_{i-1}}^{x_i} (F_x - S) = F_i - F_{i-1} - \Delta x (S_{i-1} + S_i)/2 $$

and similarly for $\Phi^{i+1/2}$. With this definition we set here

$$ \phi_i^{i-1/2} = \phi^{i-1/2}/2 + \delta A(U_{i-1/2})\phi_i^{i-1/2}, \quad \phi_i^{i+1/2} = \phi^{i+1/2}/2 - \delta A(U_{i+1/2})\phi_i^{i+1/2} \qquad (75) $$

where the term multiplied by the matrix $\delta$ is a stabilizing upwind bias, and $U_{i\pm 1/2}$ denote simple the cell averages of the unknown. Classical upwinding is obtained for $\tau = |A|^{-1}$, Lax-Wendroff like schemes are obtained for $\delta = \Delta t/2\Delta x$. This is the choice used here.

Instead of using a reconstruction, or more complex Lax-Wendroff like procedures (see e.g. [16, 11]), a second order extension is built via a two-step procedure exploiting a first update to construct a strongly consistent residual. Given the first order predictor (74) we set $(\cdot)^{n+1/2} = ((\cdot)^n + (\cdot)^{k+1})/2$, and we compute in the second stage

$$ \Phi^{i-1/2} = \frac{1}{\Delta t} \int_{t^n}^{t^{n+1}} \int_{x_{i-1}}^{x_i} (U_t + F_x - S) = \Delta x \frac{U_{i-1/2}^{k+1} - U_{i-1/2}^n}{\Delta t} + \frac{1}{2} \phi_i^{i-1/2}(U^n) + \frac{1}{2} \phi_i^{i-1/2}(U^{k+1}) $$

and similarly for $\Phi^{i+1/2}$. The update is more simply written now in the form

$$ U_i^{n+1} = U_i^{k+1} - \frac{\Delta t}{\Delta x} \Phi_i^{i-1/2} - \frac{\Delta t}{\Delta x} \Phi_i^{i+1/2} $$

where the split residuals are defined formally as in (75), replacing the $\phi$ by the $\Phi$ residuals. The interested reader can refer to [36, 1, 2] and references therein for additional information, and theoretical aspects.
Finally, in the case of the elliptic-hyperbolic approximation the time step is computed from the shallow water part as

$$ \Delta t = \nu \frac{\Delta x}{\max_i(|u_i| + \sqrt{gh_i})} \qquad (76) $$

The interested reader can refer to [11, 19, 25] for Fourier analyses justifying this choice. For the fully hyperbolic approximation the source terms are exactly integrated in the splitting method, so we use the classical condition based on the largest eigenvalue

$$ \Delta t = \nu \frac{\Delta x}{\max_i(|u_i| + \sqrt{gh_i + \mu\tau_i^2})}, \qquad (77) $$

which shows a dependence on the value of the relaxation constant $\mu$.

# 5 Numerical tests

## 5.1 Verification using analytical solutions

### 5.1.1 Solitary waves

We have performed a verification of the one-dimensional solver using the solitary wave solution of section 3. We discuss here the case of a solitary wave of nonlinearity $\epsilon = a\tau_\infty = 0.2$, and we consider the case $\chi = 0.4$ (cf. left picture on figure 4), and $h_\infty = 1m$. The behaviour with other values of the parameters is very similar.

On the spatial domain $[-200, 200]m$ we have initialized the solution centering the solitary wave at $x_0 = -100m$, and let it travel for a distance $L = 200m$. We have solved the problem both using the elliptic-hyperbolic formulation, and the hyperbolic relaxation form, and performed a grid convergence study. The results are summarized on figure 6, for both the elliptic-hyperbolic formulation (left), and for the hyperbolic one (right). Concerning the first approach, the results of the figure (and all our results) confirm a nice second order of convergence in all the norms for depth and velocity.

For the hyperbolic reformulation, the results show a dependence on the value of the relaxation constant $\mu$. Already for $\mu = 1000$ the error convergence rates reduce on finer meshes. This is due essentially to the modelling error for finite $\mu$. As the right picture on figure 6 shows, values of $\mu$ of $10^4$ are necessary to correct this. However, such values affect negatively the computational time, due to the time step restriction (77).

To clarify the impact of this on the qualitative behaviour of the model, we have looked at the convergence of the solution with $\mu$. The results are summarized in figure 7. The left plot in the figure shows the wave profile of the numerically obtained solitary wave at final time. Despite what the error convergence suggests, clearly there is no visual difference in the solutions when using values of $\mu$ above $10^3$. This is reassuring concerning the use of the hyperbolic approximation in physical applications, although the specific optimal value of $\mu$ may depend on the problem parameters. The right plot in the figure depicts the convergence of the depth wrt $\mu$, showing the expected first order rate.

### 5.1.2 Composite solution

We now consider the verification against the composite travelling solution of section §3.4. Besides verifying further the implementation, the numerical study tries to clarify the role of the regularization. As a first exercice, we perform a grid convergence study. We start from the initial solution with values provided at the end of §3.4 (cf figure 5), using the regularization with half gGN solitary wave. Simulations are run on the computational domain $[-100, 2300]m$, with the initial jump at the origin. Free boundary conditions are imposed on the right, while the constant solution is imposed on the left. We let the model evolve the flow for 50s, and then compare with the initial solution translated by the exact celerity computed analytically. The results are reported on figure 8.

To give a more general view, we report in the figure both the solutions in the subdomain $[0,575]m$, and zoom of the first peaks, including the regularized one. As one can see, both numerical approximations of the model converge nicely to the travelling solution. The difference between the original and hyperbolic reformulation are well explained by the finite value of $\mu$. Quite interestingly, the results converge to the translated initial solution, including the regularization. This suggests that the actual solution admitted is composed from a constant state, obtained from the Rankine-Hugoniot relations, and an infinite series of oscillations between the values $\tau_1$ and $\tau_2$, much resembling an infinite series of solitary waves propagating at the level $\tau_2$.

To get more insight into this we have performed the same simulation with the different regularization strategies discussed in section §3.4. Grid converged results are reported on figure 9. Clearly, past the regularized peak, the three
<figure><img src="image_8.png" /><figcaption>Figure 6: Solitary wave with $a\tau_\infty = 0.2$ and $\chi = 0.4m^4$: grid convergence. Left: elliptic-hyperbolic method. Right: hyperbolic reformulation with different values of $\mu$. Top: water height $h$. Bottom: averaged velocity.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Solitary wave with $a\tau_\infty = 0.2$ and $\chi = 0.4m^4$: convergence with $\mu$ for the hyperbolic reformulation. Left: depth distribution at final time. Right: error convergence.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: Composite solution. Grid convergence toward the regularized solution (using the gGN half soliton) at $t = 50s$ for the initial solution of figure 5 with $\chi = 0.4m^4$, $\tau_1 = 1m^{-1}$, $\tau_2 = 1.3m^{-1}$, and $\tau_3 = 1.301m^{-1}$. Top: simplified gGN model with elliptic-hyperbolic approximation. Bottom: simplified gGN in hyperbolic form with $\mu = 1000$. Left: overview of the solution in the subdomain [0,575]m. Right: first three peaks (including regularized jump).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: Composite solution: influence of the regularization of the initial solution.</figcaption></figure>

cases provide an identical infinite train of waves. Secondly, the initialization with a half solitary wave is the only case which does not introduce any secondary waves, to which the numerical approximation converges as the mesh is refined. This indicates that the exact solution is not a discontinuous one in a classical sense, but rather a composite solution composed of a constant state $\star$, and a periodic wave train of waves of large length. This traveling wave solution verifies the generalized Rankine-Hugoniot conditions coming from the mass and momentum equations. In this sense, it is a “discontinuous” solution of the gGN equations. The relation of such “discontinuous” solutions with shock solutions of the corresponding Whitham system can be established (see such a study in [22] for the Benjamin-Bona-Mahony equation having the same dispersive properties as the gGN system).

## 5.2 Breakdown of a Gaussian water column

We pass now to the study of the behaviour of the model. To begin with, we consider a simple case consisting in the breakdown of a Gaussian water column for which the initial state is given by

$$h(x, t = 0) = h_{\infty}(1 + \epsilon e^{-(x-x_0)^2/L^2}), \quad u(x, t = 0) = 0.$$

We start by investigating the propagation properties of the simplified geometrical Green-Naghdi model with respect to the dispersion and nonlinearity parameters $\chi$ and $\epsilon$. To this end we consider a battery of six cases corresponding to two values $\chi = 0.1m^4$ and $\chi = 0.6m^4$, and three non-linearity values: $\epsilon = 0.125$, $\epsilon = 0.25$, and $\epsilon = 0.5$. We consider a computational domain [-150, 150]m, on which we set the initial Gaussian centered in $x_0 = 0$, and we take $L = 2m$,
and $h_\infty = 2m$ as well. We simulate with both solution methods the evolution of the initial column up to time 16s using $N = 10^4$ cells, and taking a value $\mu = 10^3$ in the hyperbolic relaxation formulation.

The results are visualized on figure 10. Several remarks can be made. First of all, the two different formulation provide almost identical results, which is a further verification (albeit indirect) of our implementation. Second, comparing the left and right pictures, we can see that higher values of $\chi$ result in longer waves, as one may expect since in this case the dispersive regularization has more weight. The lower value of $\chi$ results in shorter and taller waves, with stronger and more localized oscillations. These differences are enhanced when considering increasing values of non-linearity. Note that as the discussion of section §2.4 shows, for a given section shape, higher values of $\chi$ may correspond both to deeper channels, as well as to larger ones.

### 5.2.1 Validation against section averaged 2D shallow water simulations

Since for this case there is no exact solution, we consider a validation against section averaged full 2D shallow water computations. The latter have been run with the solver developed in the series of works [35, 3, 4], and thoroughly validated on a wide range of cases against analytical and experimental data.

The setup of the 2D simulations is the following. We consider three sections: symmetric trapezoidal; symmetric triangular; non-symmetric trapezium. Following the notation of section §2, we set for the case $\chi = 0.1m^4 b_0 = 2.5m$, and we then for each section we choose appropriate values of the remaining lengths. In particular we take $(l_1, l_2, l_3) = (1.07, 0.36, 1.07)m$ in the symmetric trapezoidal case, and $(l_1, l_2, l_3) = (1.38575, 0.0, 1.38575)m$ in the triangular one. For the non-symmetric trapezium we choose $(l_1, l_2, l_3) = (1.284, 0.36, 0.7249)m$. For $\chi = 0.6m^4$ we have set $b_0 = 2.75m$, and chosen $(l_1, l_2, l_3) = (2.165, 1.1175, 2.165)m$ in the trapezoidal case, and $(l_1, l_2, l_3) = (3.0855, 0.0, 3.0855)m$ in the triangular one. Also in this case we have computed solutions with a non-symmetric trapezoidal section with $(l_1, l_2, l_3) = (2.598, 1.1175, 1.14668)m$.

In all cases, the initial solution is then set using equation (7):

$$h(x, y, t = 0) = \bar{h}(x, t = 0) - b(y) + \bar{b} = h_\infty(1 + \epsilon e^{-(x-x_0)^2/L^2}) - b(y) + \bar{b}.$$

with $b(y)$ defined as in figure 1m and $h_\infty = 2m$. Only the case $\epsilon = 0.25$ is commented for brevity. The computational domain is the rectangle $[-150, 150]m \times [-l_1 - l_2/2, l_2/2 + l_3]m$. Both velocity components are set to zero at $t = 0$, and reflective boundary conditions are set on all boundaries. With the above initialization, and the chosen values for $h_\infty$ and $b_0$, we satisfy all the hypotheses of the model. Fine triangulations with mesh sizes $\Delta x \approx 6.25cm$ and $\Delta y \approx 2cm$ are used. As a first step, we run a case with $b(y) = 0$. The result is reported for completeness on figure 11 showing the expected behaviour from the shallow water system with shock formation in finite time. The right plots in the figure show $y$-averaged data confirming this behaviour.

We now consider the simulations including bathymetric variations. A close up three dimensional view of the right-going waves is reported on figures 12 for the case $\chi = 0.1m^4$, and 13 for $\chi = 0.6m^4$. The figures clearly show the formation and propagation of smooth waves, with roughly a first wave resembling a solitary wave with variable amplitude along $y$, followed by secondary smooth waves. The first wave amplitudes are higher on the shallower parts than in the center. These results are similar to those presented in [12] for Favre waves: no discontinuities are produced despite of the the nonlinear hyperbolic character of the system, and of the non-linearity of the problem. Comparing the two values of $\chi$ we see that for $\chi = 0.1m^4$ (thinner or shallower channels) we obtain shorter waves, with secondary ones clustered close to the leading wave. For higher $\chi$ we see more dispersion, longer waves, and secondary waves with lower celerity. This behaviour is identical to the one of the gGN system in one dimension.

To go further, we now compute numerically $y$-averages of the 2D data, and compare them to the solutions obtained with the geometrical Green-Naghdi equations developed in this paper, which for the value of $\epsilon = 0.25$ corresponds to the second row of figure 10. The comparisons are reported in figure 14. The plots show that the averaged data for the different sections are very close and present an excellent match with the solutions of the 1D gGN model. A somewhat larger difference is observed for the non-symmetric case which clearly presents a more three dimensional wave field, with more secondary waves perhaps missed by our approximate model. The agreement is generally in favour of the new model. These results further confirm the existence of these dispersive-like waves due to transverse refraction. These are indeed dispersive waves, which can be modelled with the section averaged geometrical Green-Naghdi equations proposed in this paper.
<figure><img src="image_12.png" /><figcaption>Figure 10: Breakdown of a Gaussian water column: solutions at time $t = 16s$ for three values of the non-linearity coefficient $\epsilon = a\tau_{\infty} \in \{0.125, 0.25, 0.5\}$. Left column: $\chi = 0.1m^4$. Right column $\chi = 0.6m^4$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: Breakdown of a Gaussian water column. Free surface from shallow water simulations corresponding to $\epsilon = 0.25$. Initial solution (top-left), and solution at time $t = 5s$ (top-right). Bottom: $y$-averaged data at times $t = 5s$ (left) and $t = 15s$ (right).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 12: Breakdown of a Gaussian water column: enhanced view of the right-going wave at time $t = 8s$. Shallow water simulations on bathymetric sections corresponding to $\chi = 0.1m^4$ and $\epsilon = 0.25$. Trapezoidal (left), triangular (center), and asymmetric trapezoidal (right) sections.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: Breakdown of a Gaussian water column: enhanced view of the right-going wave at time $t = 8s$. Shallow water simulations on bathymetric sections corresponding to $\chi = 0.6m^4$ and $\epsilon = 0.25$. Trapezoidal (left), triangular (center), and asymmetric trapezoidal (right) sections.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: Breakdown of a Gaussian water column. Comparison of width averaged solutions from shallow water simulations with solutions of the gGN model. Left: $\chi = 0.1m^4$. Right: $\chi = 0.6m^4$. Solutions at time $t = 8s$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 15: Riemann problem. Comparison of solutions at $t = 15s$ for the gGN solutions using the elliptic-hyperbolic and hyperbolic relaxation for $\epsilon = 0.125$ (left column), $\epsilon = 0.25$ (central column), $\epsilon = 0.5$ (right column), and with $\chi = 0.1m^4$ (top), and $\chi = 0.6m^4$ (bottom).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 16: Riemann problem. Zoom of the leading waves for $\epsilon = 0.5$.</figcaption></figure>

## 5.3 Riemann problem

We consider now the computation of Riemann problems. As in [39] and [34] we consider a smoothed initial discontinuity:

$$h((x, t = 0)) = h_0^{\text{ID}}(x) = h_0\left(1 + \frac{\epsilon}{2}(1 - \tanh x)\right), \quad u(x, t = 0) = 0.$$

As done for the column fall, we will start by comparing the two reformulations of the model for different values of the parameters. We will consider here three values of the non-linearity $\epsilon = 0.125$, $\epsilon = 0.25$, and $\epsilon = 0.5$, and again the values of the dispersion coefficient $\chi = 0.1m^4$, and $\chi = 0.6m^4$. We compute mesh converged numerical solutions for both approaches discussed in section 4, using a value $\mu = 10^3$ for the hyperbolized formulation. The results are summarized on figure 15. The plots show that the two formulations give almost identical results, except for the highest values of the non-linarity, for which the hyperbolized model provide slightly lower peaks and a small phase shift. This effect, visible on figure 16, is most certainly due to the finite value of $\mu = 10^3$ chosen.
Concerning the qualitative behaviour of the model, as expected, smaller values of $\chi$ provide higher frequency dispersive shocks. Conversely, higher nonlinearities lead to stronger secondary waves, all travelling faster, and ultimately also dispersive shocks with shorter wavelengths at this final time.

### 5.3.1 Validation against section averaged 2D shallow water simulations

As done for the previous case, as a means of validation we have run full 2D shallow water computations. We proceed as before. We consider equivalent sections providing values of the geometrical dispersion coefficient $\chi = 0.1m^4$, and $\chi = 0.6m^4$. The geometrical definitions are the same of section §5.2.1. Then we initialize the solution based on relation (7):

$$h(x, y, t = 0) = h_0^{1D}(x) - b(y) + \bar{b},$$

with $h_0^{1D}(x)$ defined above. As before, with these values we can initialize the flow on a constant upstream average depth of 2m, without incurring in dry areas. We will consider here the case with lowest non-linearity $\epsilon = 0.125$.

To begin with, we look at the free surface obtained at $t = 15s$, which are reported on figure 17. Flat bathymetry results are reported in the bottom picture for completeness. The first thing we can see is that none of the computations with bathymetry exhibits the formation of classical shock waves, as in the flat bathymetry case. All the results show a free surface level which has an almost one dimensional phase, with finite wavelength oscillations. These wavelengths are several orders of magnitude larger than the mesh size, and have nothing to do with numerical dispersion, or lack of monotonicity, but are related to dispersion induced by physical refraction in the $y$ direction. We also note that for a given $\chi$ different sections provide very close results to one another, with perhaps the main difference between the symmetric and non-symmetric shapes. The latter ones give a more three dimensional field, richer in secondary waves. As for the water column case, the wave heights obtained are more important in shallower regions. Also here we can clearly observe, as in one dimensions, more dispersive behaviour for larger $\chi$: the waves are longer, and a larger disparity in celerity is observed.

To validate the gGN model, we compare its results with transverse averages of the 2D solutions. The resulting comparison is reported on figure 18. All the comparisons show a nice agreement. For both values of $\chi$ the symmetric 2D results provide a solution with 4 or 5 peaks, while the 1D model predicts in general more secondary oscillations. The amplitude and phases of the leading front are excellent, and those of the first two secondary waves are also captured nicely by the gGN model. The average non-symmetric results show in general slightly longer waves, which may originate from the merging of more secondary waves in this short time scale simulation.

Overall these results, as those of the previous section, validate our model, showing a very satisfactory agreement. We also confirm the findings of [12] on more simulations, and highlight the existence of these little known dispersive-like waves due to transverse refraction, which can be modelled with our geometrical Green-Naghdi equations.

### 5.3.2 Solitary wave fission

To conclude the study of the Riemann problem, we have run long time 1D simulations for the largest values $\chi = 0.6m^4$ and $\epsilon = 0.5$. The final time is now set to 600s, in order to verify the occurrence of solitary wave fission, and compare the resulting solitary wave with the ones studied in sections §3.3 and §5.1. We report in figure 19 a visualization of the front of the wave superposing several instances of the leading waves. In the visualizations, we align the second peak of the dispersive front, so that we may visualize the separation of the first solitary wave. For the last instance, we have measured an amplitude $a = 0.56h_\infty$ which we used to generate a corresponding exact gGN solitary wave solving (46) with matching amplitude. The exact solitary wave corresponding to the measured amplitude is reported in green in the picture, such that it maximum is superposed it to the first solitary wave resulting from the Riemann problem. The agreement is excellent.

## 5.4 Treske's experiments in a trapezoidal channel

As a final test, we consider the simulation of the experiments in channels with trapezoidal sections by [41]. We are interested in the issue of generation of long waves associated with geometrical dispersion. The setup of the problem is very similar to the one described in [12, 24]. The channel considered has a symmetric trapezoidal shape. Using the notation of section §3.3 (see leftmost picture on figure 2), we consider as in the experiments a width $w = 1.24m$, and a slope with $\beta = 1/3$. The initial state consists of a smoothed depth and speed discontinuity

$$\bar{h}(x) = \bar{h}_1 + \frac{\bar{h}_2 - \bar{h}_1}{2} (1 - \tanh(x/\alpha))$$

$$u(x) = \frac{u_2}{2} (1 - \tanh(x/\alpha))$$
<figure><img src="image_19.png" /><figcaption>Figure 17: Riemann problem. Free surface from shallow water solutions at $t = 15s$ for $\epsilon = 0.125$ with trapezoidal (left column) and triangular (right column) sections corresponding to $\chi = 0.1m^4$ (first row), and $\chi = 0.6m^4$ (second row). The last row shows the result for flat bathymetry.</figcaption></figure>

where, for a given Froude number, the values ($\bar{h}_2, u_2$) are obtained from section averaged Rankine-Hugoniot conditions for a bore moving on a state at rest. The value $\bar{h}_1$ is computed using the experimental height in the channel center $h_1 = 0.16m$, using the formula valid for a trapezium

$$ \bar{h} = \frac{\beta w + h}{\beta w + 2h} h. $$

To evaluate the dispersion coefficient, we use (36) with $b_0 = h_2$, the value of the central height after the bore, and with $l_2 = w$ and $l_1 = l_3 = b_0/\beta$.

The simulations are run on a domain of length $L = 100m$ using the ellipitc-hyperbolic reformation. The results of the figures are obtained on a largely converged resolution of 20.000 mesh points. As in previous works, the simulations are stopped after the first peak has travelled a distance of 72m. Finally, to choose the parameter $\alpha$ in the initialization, we use the steepness of the waves. In particular, we have set $\alpha = L(h_2 - h_1)/\lambda$ with the wavelength $\lambda$ obtained from the Lemoine theory (cf. section 3.3). We simulate the waves for values of Froude $F \in [1.0125, 1.20]$.

The left picture on figure 20 shows the typical wave profiles obtained for a low, intermediate and high value of the Froude number. We can clearly see the waves becoming shorter and steeper. As it is customary for this experiment we measure the wavelengths and amplitudes defined respectively as the distance between the first two peaks, and the height difference between the first crest and the first trough (cf. right picture on figure 20). These measures are made non-dimensional by dividing with $\bar{h}_1$, and compared to the data by [41] in figure 21. The results show a very reasonable agreement with the data and fully validate our derivation. An overprediction of the wavelength is observed for higher values of Froude. This may be related to the lack of higher order non-linear effects (and possibly dispersion) in the model. Concerning the amplitudes, as proposed in [24] we have averaged the experimental data on the channel banks and axis to obtain a fair comparison. The average is obtained as $A_{\text{banks}}/3 + 2A_{\text{axis}}/3$ corresponding to integrating a symmetric parabola. We obtain a correct prediction of the amplitude growth as well.
<figure><img src="image_20.png" /><figcaption>Figure 18: Riemann problem. Comparison of averaged free surface profiles at $t = 15s$ for $\epsilon = 0.125$, and $\chi = 0.1m^4$ (first row), and $\chi = 0.6m^4$ (second row). Geometrical Green-Naghdi model compared to transversely averaged 2D shallow water computations. Left: overall view. Right: zoom of the leading waves.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 19: Riemann problem: solitary wave fission for the case $\epsilon = 0.5$, and $\chi = 0.6m^4$. Solutions at times $t = 20s$ (blue), $t = 60s$ (black), and $t = 600s$ (red). All secondary peaks are aligned at the origin. The last solution is used to generate a corresponding exact gGN solitary wave (in green) obtained solving (46) with matching amplitude.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 20: Treske experiments in a trapezoidal channel. Left: section averaged heights obtained with the simplified gGN model for Froude numbers 1.05 (solid black), 1.1 (blue dash-dot), and 1.2 (red dashed). Right: definition of measured wavelength and amplitude.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 21: Treske experiments in a trapezoidal channel. Left: wave length as function of the Froude number. Right: amplitude as a function of the Froude number. Predictions of the simplified gGN model (squares) compared with the data by [41] (asterisks).</figcaption></figure>

## 6 Conclusion

From the two-dimensional hyperbolic non-linear shallow water equations we have derived a one-dimensional dispersive system of PDEs (geometrical Green-Naghdi type (gGN) model) describing the evolution of section averaged depth and velocity. The dispersive effects arise from transverse refraction solely due to bathymetric variations. The model proposed has several appealing properties: it is Galilean invariant; it admits an energy conservation law, and exact travelling waves solutions. Its dispersion relation allows to obtain quantitative predictions of the wavelengths of undular bores for low Froude numbers.

We have proposed two approaches to solve the model numerically, and discussed numerical tests going from comparisons with analytical travelling wave solutions, to validation with respect to multidimensional shallow water simulations, and experimental data. In particular, the multi-dimensional shallow water simulations confirm the occurrence of these dispersive-like waves in channels of finite width adding evidence to the one already presented e.g. in [12]. The section averages of the multi dimensional computations have been used to validate the geometrical gGN model proposed. Finally, the new model predicts with satisfactory accuracy the data from [41], fully confirming its ability to describe the non-linear geometrical waves observed experimentally.

Many questions remain open concerning the interaction of the different processes at play, and in particular the emergence of geometrical dispersive effects. There is certainly some the interaction between non-linearity and dispersion which controls this process, and probably some equivalent of a wave breaking range in which the shallow water equations do predict shock formation. A major issue, focus of ongoing work, is also the interaction of different dispersive
processes and the formulation of a model including with clear hypotheses both horizontal dispersion (dispersive-like waves), and usual non-hydrostatic effects as those modelled by the classical Serre-Green-Naghdi model.

Acknowledgements

SG would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme "Emergent phenomena in nonlinear dispersive waves", where work on this paper was undertaken. This work was supported by EPSRC grant EP/R014604/1. MR is a member of the CARDAMOM research team of the Inria center at University of Bordeaux.

A Geometrical coefficient properties for symmetric channels

We prove that for $b(y) = b(l - y)$ for any $0 < y < l$, i.e. for $b(y)$ symmetric with respect to the center the channel, we have

$\bar{S} = 0.$

Indeed,

$$l \bar{S} = \int_{0}^{l} \int_{0}^{y} (\bar{b} - b(s)) ds dy = \frac{\bar{b} l^2}{2} - \int_{0}^{l} \int_{0}^{y} b(s) ds dy$$

The last term is just the integration over the isosceles right triangle $A = \{(s, y) | 0 < s < y < l\}$ in the square $B = \{(s, y) | 0 < s < 1, 0 < y < l\}$. If the function $b(y)$ is symmetric, one can replace the integration over $A$ by the integration over the whole square $B = \{(s, y) | 0 < s < 1, 0 < y < l\}$ divided by 2:

$$l \bar{S} = \frac{\bar{b} l^2}{2} - \frac{1}{2} \left( \int_0^l \int_0^y b(s) ds dy + \int_0^l \int_y^l b(s) ds dy \right) = \frac{\bar{b} l^2}{2} - \frac{1}{2} \int_0^l \int_0^l b(s) ds dy \\ = \frac{\bar{b} l^2}{2} - \frac{\bar{b} l}{2} \int_0^l dy = 0.$$

B Potential for symmetric channels

We prove here proposition 1. To this end we consider the direct computation of

$$- \left( \frac{\partial N}{\partial \tau} - \frac{D}{Dt} \left( \frac{\partial N}{\partial \dot{\tau}} \right) \right) \frac{d\sigma}{dy} = \frac{d\sigma}{ds} \int_0^y \left( \frac{\sigma(s) \ddot{\tau}}{1 - \tau \frac{d\sigma(s)}{ds}} + \frac{\sigma(s) \frac{d\sigma(s)}{ds} \dot{\tau}^2}{2 \left(1 - \tau \frac{d\sigma(s)}{ds}\right)^2} \right) ds, \\ 
\frac{\partial}{\partial y} \left( \sigma(y) \frac{\partial N}{\partial \tau} \right) = \frac{d\sigma}{dy} \frac{\partial N}{\partial \tau} + \sigma(y) \frac{\partial}{\partial \tau} \left( \frac{\partial N}{\partial y} \right) \\ 
= \frac{d\sigma}{dy} \frac{\dot{\tau}^2}{2} \int_0^y \frac{\sigma(s) \frac{d\sigma(s)}{ds}}{\left(1 - \tau \frac{d\sigma}{ds}\right)^2} ds + \frac{d\sigma}{dy} \frac{\dot{\tau}^2}{2} \frac{\sigma^2(y)}{\left(1 - \tau \frac{d\sigma}{dy}\right)^2}.$$

Summing the last expressions, we obtain

$$M \frac{d\sigma}{dy} = - \left( \frac{\partial N}{\partial \tau} - \frac{D}{Dt} \left( \frac{\partial N}{\partial \dot{\tau}} \right) \right) \frac{d\sigma}{dy} + \frac{\partial}{\partial y} \left( \sigma(y) \frac{\partial N}{\partial \tau} \right). \qquad (78)$$

Now, consider the case $\bar{S} = \mathcal{O}(\varepsilon^\beta)$. This property reduces to $\bar{S} = 0$ for symmetric bathymetries $b(y)$. In this case by virtue of (15) we have

$\sigma(0) = \sigma(l) = \mathcal{O}(\varepsilon^\beta).$

Using this fact we can write

$$\overline{M \frac{d\sigma}{dy}} = - \left( \frac{\partial N}{\partial \tau} - \frac{D}{Dt} \left( \frac{\partial N}{\partial \dot{\tau}} \right) \right) \frac{d\sigma}{dy} + \sigma(y) \frac{\partial N}{\partial \tau} \Big|_{0}^{l} \qquad (79)$$
$$ = -\left(\frac{\partial N}{\partial \tau} - \frac{D}{Dt} \left(\frac{\partial N}{\partial \dot{\tau}}\right)\right) \frac{d\sigma}{dy} + \mathcal{O}(\varepsilon^\beta) = -\left(\frac{\partial N \frac{d\sigma}{dy}}{\partial \tau} - \frac{D}{Dt} \left(\frac{\partial N \frac{d\sigma}{dy}}{\partial \dot{\tau}}\right)\right) + \mathcal{O}(\varepsilon^\beta). \quad (80) $$

Denoting

$$ \mathcal{L}(\tau, \dot{\tau}) = N \frac{d\sigma}{dy} $$

we finally get in the case of quasi-symmetric channels, and within $\mathcal{O}(\varepsilon^\beta)$:

$$ \overline{M \frac{d\sigma}{dy}} = - \left( \frac{\partial \mathcal{L}}{\partial \tau} - \frac{D}{Dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{\tau}} \right) \right) = - \frac{\delta \mathcal{L}}{\delta \tau} \quad (81) $$

## References

[1] R. Abgrall and M. Ricchiuto. High order methods for CFD. In René de Borst Erwin Stein and Thomas J.R. Hughes, editors, *Encyclopedia of Computational Mechanics*, Second Edition. John Wiley and Sons, 2017.

[2] Rémi Abgrall and Mario Ricchiuto. *Hyperbolic Balance Laws: Residual Distribution, Local and Global Fluxes*, pages 177–222. Springer Nature Singapore, Singapore, 2022.

[3] Luca Arpaia and Mario Ricchiuto. r- adaptation for shallow water flows: conservation, well balancedness, efficiency. *Computers & Fluids*, 160:175–203, 2018.

[4] Luca Arpaia and Mario Ricchiuto. Well balanced residual distribution for the ALE spherical shallow water equations on moving adaptive meshes. *J.Comput.Phys.*, 405, 2020.

[5] A. Berezovski, J. Engelbrecht, and M. Berezovski. Waves in microstructured solids: a unified viewpoint of modeling. *Acta Mech.*, 220:349–363, 2011.

[6] A. Berezovski, J. Engelbrecht, A. Salupere, K. Tamm, T. Peets, and M. Berezovski. Dispersive waves in microstructured solids. *International Journal of Solids and Structures*, 50(11):1981–1990, 2013.

[7] N. Bonneton, P. Bonneton, J.-P. Parisot, A Sottolichio, and G Detandt. Tidal bore and Mascaret - example of Garonne and Seine Rivers. *Comptes Rendus Geosciences*, 344(508-515), 2012.

[8] P. Bonneton, N. Bonneton, J.-P. Parisot, and B. Castelle. Tidal bore dynamics in funnel-shaped estuaries. *Journal of Geophysical Research: Oceans*, 120(2):923–941, 2015.

[9] Philippe Bonneton, Florent Chazel, David Lannes, Fabien Marche, and Marion Tissier. A splitting approach for the fully nonlinear and weakly dispersive Green–Naghdi model. *J.Comput.Phys*, 230(4):1479–1498, 2011.

[10] S. Busto, M. Dumbser, C. Escalante, N. Favrie, and S. Gavrilyuk. On High Order ADER Discontinuous Galerkin Schemes for First Order Hyperbolic Reformulations of Nonlinear Dispersive Systems. *Journal of Scientific Computing*, 87, 2021.

[11] Aurore Cauquis, Mario Ricchiuto, and Philippe Heinrich. Lax–wendroff schemes with polynomial extrapolation and simplified lax–wendroff schemes for dispersive waves: A comparative study. *Water Waves*, pages 1–33, 2022.

[12] R. Chassagne, A.G. Filippini, M. Ricchiuto, and P. Bonneton. Dispersive and dispersive-like bores in channels with sloping banks. *Journal of Fluid Mechanics*, 870:595–616, 2019.

[13] Rémi Chassagne, Andrea Gilberto Filippini, Mario Ricchiuto, and Philippe Bonneton. Dispersive and dispersive-like bores in channels with sloping banks. Research Report RR-9228, Inria Bordeaux Sud-Ouest, 2018.

[14] Firas Dhaouadi, Nicolas Favrie, and Sergey Gavrilyuk. Extended Lagrangian approach for the defocusing nonlinear Schrödinger equation. *Studies in Applied Mathematics*, 142(3):336–358, 2019.

[15] Vincent Duchêne. Rigorous justification of the Favrie–Gavrilyuk approximation to the Serre–Green–Naghdi model. *Nonlinearity*, 32(10):3772, sep 2019.

[16] Michael Dumbser, Dinshaw S. Balsara, Eleuterio F. Toro, and Claus-Dieter Munz. A unified framework for the construction of one-step finite volume and discontinuous galerkin schemes on unstructured meshes. *J.Comput.Phys.*, 227(18):8209–8253, 2008.

[17] H. Favre. *Etude théorique et expérimentale des ondes de translation dans les canaux découverts*. Dunod, 1935.
[18] N Favrie and S Gavrilyuk. A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves. *Nonlinearity*, 30(7):2718, 2017.

[19] Andrea Gilberto Filippini, Maria Kazolea, and Mario Ricchiuto. A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up. *J.Comput.Phys.*, 310:381–417, 2016.

[20] Sergey Gavrilyuk. Multiphase flow modeling via Hamilton’s principle. In Francesco dell’Isola and Sergey Gavrilyuk, editors, *Variational Models and Methods in Solid and Fluid Mechanics*, pages 163–210. Springer Vienna, Vienna, 2011.

[21] Sergey Gavrilyuk, Boniface Nkonga, Keh-Ming Shyue, and Lev Truskinovsky. Stationary shock-like transition fronts in dispersive systems. *Nonlinearity*, 33(10):5477, 2020.

[22] Sergey Gavrilyuk and Keh-Ming Shyue. Singular solutions of the BBM equation: analytical and numerical study. *Nonlinearity*, 35(1):388, nov 2021.

[23] Sergey Gavrilyuk and Keh-Ming Shyue. 2D Serre-Green-Naghdi equations over topography: Elliptic operator inversion method. *Journal of Hydraulic Engineering*, 150(1), 2024.

[24] B. Jouy, D. Violeau, M. Ricchiuto, and M. Le. One dimensional modelling of Favre waves in channels. *Applied Mathematical Modelling*, 133:170–194, 2024.

[25] M. Kazolea, A.G. Filippini, and M. Ricchiuto. Low dispersion finite volume/element discretization of the enhanced green-naghdi equations for wave propagation, breaking and runup on unstructured meshes. *Ocean Modelling*, 182:102157, 2023.

[26] D. Ketcheson and M. Quezada de Luna. Diffractons : Solitary Waves Created by Diffraction in Periodic Media. *Multiscale Modeling & Simulation*, 13(1):440–458, January 2015.

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

[28] R. Lemoine. Notules hydrauliques. *La Houille Blanche*, (2):183–186, March 1948.

[29] Alireza Mazaheri, Mario Ricchiuto, and Hiroaki Nishikawa. A first-order hyperbolic system approach for dispersion. *J.Comput.Phys.*, 321:593–605, 2016.

[30] Victor Michel-Dansac, Pascal Noble, and Jean-Paul Vila. Consistent section-averaged shallow water equations with bottom friction. *Eur.J.Mech. - B/Fluids*, 86, 2021.

[31] M. Parisot. Entropy-satisfying scheme for a hierarchy of dispersive reduced models of free surface flow. *International Journal for Numerical Methods in Fluids*, 91(10):509–531, 2019.

[32] D.H. Peregrine. Long waves in a uniform channel of arbitrary cross-section. *Journal of Fluid Mechanics*, 32:353–365, 1968.

[33] Manuel Quezada de Luna and David I. Ketcheson. Solitary water waves created by variations in bathymetry. *Journal of Fluid Mechanics*, 917, May 2021.

[34] Hendrik Ranocha and Mario Ricchiuto. Structure-preserving approximations of the Serre-Green-Naghdi equations in standard and hyperbolic form, 2024.

[35] M. Ricchiuto. An explicit residual based approach for shallow water flows. *J.Comput.Phys.*, 280:306–344, 2015.

[36] Mario Ricchiuto. *Contributions to the development of residual discretizations for hyperbolic conservation laws with application to shallow water flows*. Habilitation à diriger des recherches, Université Sciences et Technologies - Bordeaux I, December 2011.

[37] M. H. Teng and T. Y. Wu. Nonlinear water waves in channels of arbitrary shape. *Journal of Fluid Mechanics*, 242:211–233, 1992.

[38] M. H. Teng and T. Y. Wu. Effects of channel cross-sectional geometry on long wave generation and propagation. *Physics of Fluids*, 9:3368–3377, 1997.

[39] Sergey Tkachenko, Sergey Gavrilyuk, and Jacques Massoni. Extended Lagrangian approach for the numerical study of multidimensional dispersive waves: Applications to the Serre-Green-Naghdi equations. *J.Comput.Phys.*, 477:111901, 2023.

[40] Davide Torlo and Mario Ricchiuto. Model order reduction strategies for weakly dispersive waves. *Mathematics and Computers in Simulation*, 205:997–1028, 2023.

[41] A. Treske. Undular bores (Favre-waves) in open channels - Experimental studies. *Journal of Hydraulic Research*, 32(3):355–370, May 1994.
[42] Ge Wei, James T Kirby, Stephan T Grilli, and Ravishankar Subramanya. A fully nonlinear Boussinesq model for surface waves. part 1. highly nonlinear unsteady waves. *Journal of Fluid Mechanics*, 294:71–92, 1995.

[43] Patricio Winckler and Philip L.-F. Liu. Long waves in a straight channel with non-uniform cross-section. *Journal of Fluid Mechanics*, 770:156–188, May 2015.