# One dimensional modelling of Favre waves in channels

Bastien Jouy, Damien Violeau, Mario Ricchiuto, Ming Hoang Le

► To cite this version:

Bastien Jouy, Damien Violeau, Mario Ricchiuto, Ming Hoang Le. One dimensional mod-
elling of Favre waves in channels. Applied Mathematical Modelling, 2024, 133, pp.170-194.
10.1016/j.apm.2024.05.020. hal-04706732

HAL Id: hal-04706732

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

Submitted on 23 Sep 2024

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
# One dimensional modelling of Favre waves in channels

B. Jouy$^{a,b}$, D. Violeau$^{a,b}$, M. Ricchiuto$^c$, M. Le$^b$

$^a$EDF R&D / LNHE, 6 quai Watier, Chatou, 78400, France

$^b$Saint-Venant Laboratory for Hydraulics, 6 quai Watier, Chatou, 78400, France

$^c$Inria Bordeaux - Sud-Ouest, 200 Avenue de la Vieille Tour, Talence, 33405, France

## Abstract

In this study, we propose a modified version of section-averaged Boussinesq equations of Winckler-Liu. The model is reformulated in conservative variables, allowing a decoupling of the Shallow-Water equations and the dispersive problem. An appropriate hybrid finite volume and finite element discretisations are performed and verified with a solitary wave solution derived for the typical case of prismatic channels with a trapezoidal cross-section. Classical upwind and energy conservative schemes are investigated through channel Favre waves cases. The proposed model and numerical approximations can accurately reproduce the main features of the wave train's free surface. Numerical dissipation introduced by decentered flux formulations is discussed, emphasising the need for precaution in their applications to realistic engineering issues.

**Keywords:** Favre waves, Dispersive waves, Boussinesq-type model, Finite Volume, Finite Element, Entropy Conservative

## 1. Introduction

The numerical discretization of Boussinesq-type models and more generally dispersive wave equations, including scalar models like Korteweg-de
Vries (KdV) and Benjamin-Bona-Mahony (BBM) equations, is crucial in numerically modelling wave propagation phenomena. These equations, commonly employed to describe wave characteristics in dispersive mediums such as channels or coastal areas, exhibit nonlinear and dispersive features that require careful consideration regarding discretization. Over the years, numerous methods for expressing equations to facilitate their numerical approximation have been proposed. Through these reformulations, various discretization approaches have been applied, encompassing the most conventional methods such as finite differences (FD), finite volumes (FE), finite elements (FE), and combinations of the three. The choice of a specific discretization method depends on the characteristics of the equations being solved and the physical scenarios being modelled. This involves striking a balance between accuracy, stability, and computational efficiency, ensuring that the chosen discretization scheme effectively represents the nonlinear and dispersive behaviours inherent in the modelled wave phenomena.

The seminal work by [1] employs an unsplit approach to discretize each term of the partial differential equation (PDE) system proposed by Nwogu (1993) using finite differences of a suitable order. The discrete equations obtained have time derivative terms like $G_t$, where $G = \mathcal{L}(u)$, with $u$ representing the velocity and $\mathcal{L}$ being a (usually elliptic) differential operator. This unsplit method has been the basis of numerous works extending its principles, particularly in developing finite element methods (FEM) [2–7] and finite volume [8, 9] approximations. In [10], a general setting embedding both finite differences and finite elements (continuous and discontinuous) is proposed.

Various approaches aimed at addressing the dispersive terms somewhat in-
dependently in the context of solving the Shallow Water (SW) equations. In [11], a classical time-splitting technique is employed to solve the SW equations. The work of [12] introduces a differential formulation that explicitly incorporates an algebraic source term into the SW equations. To retrieve the sought dispersive equation, the source term is shown to be constrained by an appropriately defined second-order differential problem. In [13], a discrete variant of the latter approach is performed, using the concept of projection in divergence-constrained spaces, which has been extensively developed in the past for incompressible equations. Variations of these concepts have been further explored in works such as [14], [15], [16], [17], and [18]. These subsequent studies leverage finite element, finite difference, finite volume, and hybrid approaches to discretize Boussinesq-type equations, showcasing the versatility and adaptability of these foundational ideas across different numerical methods.

Another approach consists of hyperbolic reformulations of equations suitable for applying numerical methods similar to those used to solve the non-dispersive shallow water model. In [19], a generalization of Cattaneo's original work on the heat equation is presented. This method focuses on solving the original dispersive system, using a fully hyperbolic relaxation at each time step to simultaneously approximate both the elliptic operator and the underlying non-dispersive part. The method inherently adopts an implicit time-stepping approach. A second approach introduced by [20] proproses a new augmented first-order system with stiff algebraic relaxation terms. The original dispersive model is recovered in the limit of instantaneous (infinitely stiff) relaxation. This method allows for explicit or semi-implicit time step-
ping but introduces a small modelling error. Discretization methods based
on finite volumes, finite differences, or high-order finite elements and exploit-
ing this equation form have been proposed in works such as [21], [22], [23],
and [24].

In the realm of weakly nonlinear and weakly dispersive waves within channels that are not too broad, incorporating changes in channel geometry becomes crucial for a more accurate representation of wave propagation. These types of waves are commonly encountered in engineering applications, such as tidal bores [25] or during emergency operations involving sluice gates or turbines in hydropower plants [26, 27]. In the latter scenario, the triggering induces the propagation of a dispersive shock wave known as Favre waves [28], capable of propagating several kilometres and exhibiting complex dynamics that may lead to flooding risks. The propagation dynamics are intricately governed by the shape of the channel’s cross-section [29, 30], necessitating consideration in the selected equations.

The modelling of weakly dispersive waves in channels with variable cross-sections has been the subject of many research, encompassing Boussinesq-type and KdV equations [31–38]. Despite these efforts, mentioned theories fall in adequately addressing simultaneously arbitrary cross-sections and substantial changes in channel geometry within the wavelength scale in the direction of wave propagation. Recent advancements include the development of models that overcome these limitations, employing traditional perturbation expansion methods for the Boussinesq equations [39] and also the fully nonlinear Serre equations [40]. Alternatively, a variational approach has been used, leading to a straightforward derivation of an extension of the Serre
equations tailored for variable cross-sections [41].

Our main aim is to study the numerical approximation of the section averaged Boussinesq type equations, and in particular the model by Winckler and Liu [39]. First, we propose a modified variant of the model based on the exact mass conservation equation. We then adopt the approach proposed by [12] and recast the dispersive system as the section averaged SW system plus a source, which we show being the solution of a second order (elliptic) differential operator, which we demonstrate to be coercive. We then propose a hybrid finite volume/element approximation of the resulting system. More particular, a second order $P_1$ FEM is used to approximate the elliptic operator providing the dispersive source term, and a conservative and well-balanced finite volume approximation is used for the shallow water sub-system. The numerical model is verified using a solitary solution derived from this reformulated conservative variables model. Particular attention is given to the impact of the dissipation in the FV solution equations. To this end, consider both Roe schemes of orders one, two and three, as well as an extension of the Energy Conservative (EC) scheme by [42] to the case of the section averaged SW equations developed here. Applied to the Favre waves cases, our results indicated that for long times propagation, the artificial dissipation introduced by an upwind scheme, such as Roe's formulation had to be considered carefully and treated with high-order (at least 4th) accuracy. In contrast, the EC scheme gives satisfactory results with coarse meshes. We then validate the model with laboratory experiments of [28] and [30], showing good agreements for the maximum amplitude of the leading wave. However, we also show the need to consider the temporal dependencies of the char-
acteristic wavelength to get valid comparisons with experimental data. The absence of significant effects is evident in the various numerical flux formulations, and adding friction effects does not yield any noteworthy impact at these channel scales. Therefore, caution is advised when applying numerical models at larger scales in the context of these validations.

## 2. Governing Equations

### 2.1. Generalities and notation

To model the free surface flow in a channel we consider a reference system with spatial axis $x$ aligned with the flow direction, and $(y, z)$ plane describing the generic channel section. In this system we introduce the notations sketched in figure 1. In particular, we denote by $\eta = \eta(\mathbf{x}, t)$ the elevation of the free surface wrt a reference depth $h_0(\mathbf{x})$, where $t$ denotes the time, and $\mathbf{x} = (x, y)^T$ the horizontal coordinates. Similarly, we denote by $h = h(\mathbf{x}, t)$ the water depth which verifies the identity

$$h(\mathbf{x}, t) = h_0(\mathbf{x}) + \eta(\mathbf{x}, t). \qquad (2.1)$$

Next, we define $A_0(x)$ the domain of wetted cross section at rest and we denote $A_0(x) \doteq |A_0(x)|$ the wetted area corresponding to the local depth at rest $h_0(\mathbf{x})$. We also denote by $B_0(x) \doteq |B_0(x)|$ the corresponding width of free surface at rest $B_0(x)$. The total wet area, accounting for the total depth $h(\mathbf{x}, t)$, is instead denoted by $A(x, t)$. We can now introduce the section and
width averages

$$ \begin{aligned} \langle (\cdot)(x, y, z, t) \rangle &\doteq \frac{1}{A(x, t)} \int_{A(x,t)} (\cdot) dy dz, \\ \overline{(\cdot)(x, y, z, t)} &\doteq \frac{1}{A_0(x)} \int_{A_0(x)} (\cdot) dy dz, \end{aligned} \qquad (2.2) $$

and

$$ (\cdot) \tilde{(x, y, t)} \doteq \frac{1}{B_0(x)} \int_{B_0(x)} (\cdot) dy. \qquad (2.3) $$

Both section and width averaged quantities are only functions of $x$ and $t$.

Finally, we consider reference values of the wave amplitude $\hat{a}_0$, wavelength $\hat{\lambda}$, as well as a reference (constant) depth $\hat{h}_0$. These values allow to define the dimensionless numbers

$$ \mu = \frac{\hat{h}_0}{\hat{\lambda}}, \quad \epsilon = \frac{\hat{a}_0}{\hat{h}_0}, \qquad (2.4) $$

which measure respectively the frequency dispersion ratio and the wave non-linearity.

## 2.2. Winckler-Liu's equations

We start here from the asymptotic cross-sectional model proposed by Winckler and Liu [39], which in dimensionless variables can be written as:

$$ \tilde{\eta}_t + \frac{1}{B_0} (\langle u \rangle A_0)_x - \epsilon \frac{B_0'}{B_0^2} \tilde{\eta} (\langle u \rangle A_0)_x + \epsilon \frac{1}{B_0} (\langle u \rangle B_0 \tilde{\eta})_x = O(\epsilon^2, \epsilon\mu^2, \mu^4), \quad (2.5) $$

$$ \langle u \rangle_t + \epsilon \langle u \rangle \langle u \rangle_x + g \tilde{\eta}_x + \mu^2 [\alpha \langle u \rangle_t + \beta \langle u \rangle_{tx} + \gamma \langle u \rangle_{txx}] = O(\epsilon^2, \epsilon\mu^2, \mu^4), \quad (2.6) $$

where $u(x, y, z, t)$ is the horizontal flow velocity in the $x$ direction, and $B_0'$ denotes the sidewall slope. The subscripts $x$ and $t$ denote partial derivatives
<figure><img src="image_2.png" /><figcaption>Figure 1: Definition sketch of the fluid domain: (a) in the x-z plan, and (b) in the y-z plan.</figcaption></figure>

in the $x$ direction and in time. The coefficients $\alpha(x)$, $\beta(x)$ and $\gamma(x)$ are channel geometry dependent. They are defined as

$$ \alpha = (\tilde{\chi}_1)_x - \overline{\chi_{1,x}}, \qquad (2.7) $$

$$ \beta = \tilde{\chi}_1 - \overline{\chi}_1 + (\tilde{\chi}_2)_x - \overline{\chi}_{2,x}, \qquad (2.8) $$

$$ \gamma = \tilde{\chi}_2 - \overline{\chi}_2, \qquad (2.9) $$

with $\chi_1$ and $\chi_2$ solutions of the following boundary value problems (BVPs):

$$ \Delta\chi_1 = 0, \qquad \text{over } \mathcal{A}_0, \qquad (2.10) $$

$$ \nabla\chi_1 \cdot \mathbf{n} = -\frac{A_{0,x}}{B_0}, \qquad \text{on } \mathcal{B}_0, \qquad (2.11) $$

$$ = -\frac{h_{0,x}}{\sqrt{h_{0,y}^2 + 1}}, \qquad \text{on } \partial\mathcal{A}_0\setminus\mathcal{B}_0, \qquad (2.12) $$
and

$$ \Delta\chi_2 = -1, \qquad \text{over } \mathcal{A}_0, \qquad (2.13) $$

$$ \nabla\chi_2 \cdot \mathbf{n} = -\frac{A_0}{B_0}, \qquad \text{on } \mathcal{B}_0, \qquad (2.14) $$

$$ = 0, \qquad \text{on } \partial\mathcal{A}_0\setminus\mathcal{B}_0, \qquad (2.15) $$

where $\mathbf{n}$ is a normal unit outward vector to $\partial\mathcal{A}_0$ and $\partial\mathcal{A}_0\setminus\mathcal{B}_0$ denotes the wetted perimeter (i.e. the banks and the bottom). The main geometrical notation is depicted in Figure 1.

## 2.3. Conservative reformulation of Winckler-Liu's model

Winckler-Liu's equations (2.5)-(2.6) are expressed in terms of non-conversative variables, namely $\tilde{\eta}$ and $\langle u \rangle$. In this section, we propose a reformulation allowing to work with the total wet cross-section $A$ and the total flow rate $Au$. Note that, for the sake of simplicity, we will now omit to explicit denote averaging operators in $\langle u \rangle$ and $\tilde{\eta}$, while still working with section averaged velocities, and width averaged elevation.

As a first step, we note that (2.5) is an approximation of the exact volume conservation equation (see [41] and references therein) :

$$ A_t + \epsilon (Au)_x = 0. \qquad (2.16) $$

Next, we introduce the elliptic operator $\mathcal{T}$ defined as

$$ \mathcal{T} \doteq \alpha \mathcal{I} + \beta^* \partial_x - \partial_x (\gamma^* \partial_x), \qquad (2.17) $$

where $\gamma^* \doteq -\gamma$ and $\beta^* \doteq \beta + \gamma_x^*$, with $\mathcal{I}$ the identity operator. We can now rewrite the momentum equation (2.6) as follows:

$$ (1 + \mu^2 \mathcal{T}) u_t + \epsilon u u_x + g \eta_x = O(\epsilon^2, \epsilon \mu^2, \mu^4). \qquad (2.18) $$
Multiplying the continuity equation (2.16) by $u$ and the momentum equation (2.18) by $A$, and then summing the two, the momentum equation becomes

$$ (Au)_t + \epsilon (Au^2)_x + gA\eta_x + \mu^2 A\mathcal{T}u_t = O(\epsilon^2, \epsilon\mu^2, \mu^4). \quad (2.19) $$

In order to put the previous system into conservative form, we define the new operator

$$ \mathcal{T}(\cdot) \doteq A\mathcal{T}\left(\frac{1}{A}(\cdot)\right). \quad (2.20) $$

After a few manipulations and using the continuity equation (2.16), one can easily show that the last term of the left-hand side (LHS) of (2.19) can be written as

$$ A\mathcal{T}u_t = \mathcal{T}[(Au)_t] + \mathcal{T}[\epsilon u(Au)_x]. $$

We now consider the section averaged shallow water residual $\phi$ satisfying

$$ A\phi \doteq (Au)_t + \epsilon (Au^2)_x + gA\eta_x. \quad (2.21) $$

Denoting by $H$ the hydrodynamic head

$$ H = \epsilon \frac{u^2}{2} + g\eta, \quad (2.22) $$

we observe that

$$ A(\epsilon u u_x + g\eta_x) = A\left(\epsilon \frac{u^2}{2} + g\eta\right)_x = AH_x. \quad (2.23) $$

Thus, applying the operator $1 + \mu^2 \mathcal{T}$ to the Saint-Venant residual (2.21) and subtracting the result to (2.19) we obtain

$$ (1 + \mu^2 \mathcal{T})(A\phi) = \mu^2 \mathcal{T}(AH_x). \quad (2.24) $$
Note that the term $T(AH_x)$ involves third-order derivatives. To avoid this, by making the change of variable $\Phi \doteq \phi - H_x$ we rewrite this equation in the following form, using the definition (2.20):

$$ (1 + \mu^2 T) \Phi = -H_x. \qquad (2.25) $$

On the other hand, the momentum equation can now be rewritten as:

$$ (Au)_t + \epsilon (Au^2)_x + gA\eta_x = A\Phi + AH_x. \qquad (2.26) $$

In addition, this system can be put into conservative form by expanding the term involving the derivative of the free surface elevation:

$$ (Au)_t + \epsilon (Au^2 + gK)_x = A\Phi + AH_{s,x}, \qquad (2.27) $$

where

$$ K(h) \doteq \int A(h) dh + C, \qquad (2.28) $$

with C a constant and,

$$ H_s \doteq \epsilon \frac{u^2}{2} + gh \qquad (2.29) $$

is the hydrostatic specific head. The resulting system is composed of the continuity equation (2.16), the momentum equation (2.27) and the additional elliptic equation (2.25):

$$ \left\{ \begin{array}{l} (1 + \mu^2 T) \Phi = -H_x, \\ A_t + \epsilon (Au)_x = 0, \\ (Au)_t + (\epsilon Au^2 + gK)_x = A\Phi + AH_{s,x}. \end{array} \right. \qquad (2.30) $$

This formulation, allows to write the section averaged dispersive model as a perturbation of the section averaged Saint-Venant equations, with a dispersive forcing term obtained from a time independent differential operator
which can be solved in a decoupled manner. This formulation follows the one used to numerically solve for the Serre-Green-Naghdi (SGN) and other Boussinesq models in [12, 16-18, 43]. More precisely, we will consider a resolution method in which we solve the Saint-Venant equations with the additional forcing terms involving $\Phi$. However, when solving the elliptic equation, the operator $1 + \mu^2\mathcal{T}$ must be inverted and its inverse must therefore be bounded in order to obtain uniqueness of the solution, according to the Lax-Milgram theorem. This is addressed in the next section.

## 2.4. Coercivity of the elliptic operator (prismatic channel case)

To study the coercivity of the operator $1 + \mu^2\mathcal{T}$, we look at the quadratic form $a(u, u)$ associated with the variational formulation of the elliptic equation (2.25). Multiplying by a test function $v(x)$ and integrating over $\mathbb{R}$, after integration by parts of the first term, we obtain:

$$a(u, u) \doteq \mu^2 \int_{\mathbb{R}} \gamma^* u_x^2 dx + \int_{\mathbb{R}} \left( \alpha^* - \frac{1}{2} \mu^2 \beta_x^* \right) u^2 dx, \quad (2.31)$$

where $\alpha^* \doteq 1 + \mu^2\alpha$. Hence, $\alpha^* - \frac{1}{2}\mu^2\beta_x^* > 0$ and $\gamma^* > 0$ are sufficient conditions required for $1 + \mu^2\mathcal{T}$ to be coercive. At this stage, we could not demonstrate the first inequality theoretically. In the following, we consider only prismatic channels, i.e. the undisturbed wetted cross-section $A_0$ is independent of $x$. In this case, $\alpha = \beta = 0$ and $\gamma^* = \text{cst}$ [39]. We start with the expression of $\gamma^*$ (2.9) to prove its positivity in this case:

$$\gamma^* = \frac{1}{A_0} \int_{A_0} \chi_2 \, dS - \frac{1}{B_0} \int_{B_0} \chi_2 \, d\ell, \quad (2.32)$$

where $d\ell = \sqrt{1 + h_{0,y}^2} dy$ is a length element of $\partial A_0$ and $dS = dydz$ is a surface element of $A_0$. Using Poisson's equation (2.13) and Leibniz's rule for
the first integral, we obtain:

$$ \gamma^* = -\frac{1}{A_0} \int_{A_0} \nabla \cdot (\chi \nabla \chi) \, dS + \frac{1}{A_0} \int_{A_0} |\nabla \chi|^2 \, dS - \frac{1}{B_0} \int_{B_0} \chi \, d\ell. $$

We now use the divergence theorem for the first integral and the notation corresponding to the average over $A_0$ for the second:

$$ \gamma^* = -\frac{1}{A_0} \int_{\partial A_0} \chi \nabla \chi \cdot n \, d\ell + \overline{|\nabla \chi|^2} - \frac{1}{B_0} \int_{B_0} \chi \, d\ell. $$

Finally, we use the boundary conditions (BCs) (2.14) and (2.15) for the integral over $\partial A_0$, and we arrive at the following relationship:

$$ \gamma^* = \overline{|\nabla \chi|^2} > 0. \qquad (2.33) $$

The coercivity of the elliptic operator is therefore proved for any section $A_0$ in a prismatic channel. In [39], analytical expressions of $\gamma$ are given for prismatic channels of rectangular and triangular cross-section. For the latter geometry the author found that the sign of $\gamma$ depends on a shape ratio between $B_0$ and $h_0$. Our equation (2.33) shows that this result is wrong; our amended analytical expression for $\gamma^*$ is given in Appendix A.

### 3. Numerical approximation

In this section we focus on the numerical solving of the system (2.30). As mentioned in section 2, decoupling allows us to treat the hyperbolic and elliptic parts separately. As in [12], we use a hybrid finite element-finite volume method for the elliptic and hyperbolic part, respectively.
### 3.1. Hyperbolic part - Finite Volume Method

The continuity equation (2.16) and the momentum equation (2.27) amounts to a system of balance laws, written in dimensional form as:

$$ \mathbf{U}_t + \mathbf{F}(\mathbf{U})_x = \mathbf{S}(\mathbf{U}), \qquad (3.1) $$

with,

$$ \mathbf{U} \doteq \begin{pmatrix} A \\ Au \end{pmatrix}, \quad \mathbf{F} \doteq \begin{pmatrix} Au \\ Au^2 + gK(h) \end{pmatrix} \quad \text{and} \quad \mathbf{S} \doteq \begin{pmatrix} 0 \\ A\Phi + AH_{s,x} \end{pmatrix}, \qquad (3.2) $$

where $\mathbf{U}$, $\mathbf{F}$, $\mathbf{S}$ are the vectors of the unknown conservative variables, fluxes and source terms, respectively. Denoting the cell $[x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}]$, and for simplicity a constant grid size $\Delta x \doteq x_{i+\frac{1}{2}} - x_{i-\frac{1}{2}}$, we define the cell average of $\mathbf{U}$ as follows:

$$ \mathbf{U}_i(t) \doteq \frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \mathbf{U}(x, t) \, dx. \qquad (3.3) $$

We now use the finite volume method, which is similar in the notation to that used in [12], and is written as follows:

$$ \frac{d}{dt} \mathbf{U}_i = -\frac{1}{\Delta x} \left( \mathbf{F}_{i+\frac{1}{2}} - \mathbf{F}_{i-\frac{1}{2}} \right) + \frac{1}{\Delta x} \Delta \mathbf{S}_{H_{s,i}} + \bar{\mathbf{S}}_{\Phi}, \qquad (3.4) $$

where $\mathbf{F}_{i \pm \frac{1}{2}}$, $\Delta \mathbf{S}_{H_{s,i}}$ and $\bar{\mathbf{S}}_{\Phi}$ are the numerical fluxes and source terms at cell interface $x = i \pm \frac{1}{2}$, respectively.

These numerical fluxes can be expressed in various ways, e.g. Roe, Harten-Lax-van Leer (HLL) (see e.g. the classical book [44] for a review). Differently from the hyperbolic shallow water equations, dispersive Boussinesq-type equations do not develop genuine shock waves. Indeed, even undular
bores as the well known Favre waves [27, 28, 43] are essentially smooth waves which can propagate over very long distances. It is also known that the presence of physical dissipation processes has an impact on the dynamics of such waves, especially for long time propagation [45, 46]. Unless additional modelling terms are included such processes are not accounted for in equations (3.1). However, in practical computations numerical dissipation may replace such models, and considerably impact the wave dynamics predicted when solving numerically (3.1), or any other Boussinesq-type model. To investigate this aspect, in this work we have decided to use two types of numerical fluxes: the classical upwind Roe flux [47], and the exactly energy conservative (EC) scheme by [42]. The latter is exactly energy conservative in absence of dispersive terms.

### 3.1.1. Roe scheme

The numerical flux of (3.4) can be written in Roe's formulation as:

$$F_{i+\frac{1}{2}} = \frac{1}{2} [F(U_{i+\frac{1}{2}}^L) + F(U_{i+\frac{1}{2}}^R)] - \frac{1}{2} |\tilde{J}|_{i+\frac{1}{2}} \Delta \tilde{U}_{i+\frac{1}{2}}, \quad (3.5)$$

where $\Delta \tilde{U}_{i+\frac{1}{2}} \doteq U_{i+\frac{1}{2}}^R - U_{i+\frac{1}{2}}^L$ and $U_{i+\frac{1}{2}}^L, U_{i+\frac{1}{2}}^R$ are the state vectors on the left and right sides of interface $i + \frac{1}{2}$, respectively. They are computed using 2nd and 3rd-order MUSCL reconstruction [48]. Also, $\tilde{J}_{i+\frac{1}{2}}$ denotes the Jacobian matrix of system (3.4) in the Roe average state:

$$\tilde{J}_{i+\frac{1}{2}} = (\tilde{P}\tilde{\Lambda}\tilde{P}^{-1})_{i+\frac{1}{2}}, \quad (3.6)$$

with the eigenvectors matrices $\tilde{P}$ and the diagonal matrix of eigenvalues $\tilde{\Lambda}$. Eigenvalues are expressed from average velocity $\tilde{u}$ and celerity $\tilde{c}$ as follows:

$$\tilde{\lambda}_{1,2} = \tilde{u} \pm \tilde{c}, \quad \tilde{u}_{i+\frac{1}{2}} = \frac{u^L \sqrt{A^L} + u^R \sqrt{A^R}}{\sqrt{A^L} + \sqrt{A^R}}, \quad \tilde{c}_{i+\frac{1}{2}} = \sqrt{\frac{1}{2}g(\bar{h}^L + \bar{h}^R)}, \quad (3.7)$$
with $\bar{h} \doteq \frac{A}{B}$, the mean depth water.

The discretisation of the source terms plays an important role in the construction of a numerical scheme. In particular, it has been shown that a projection of the source terms onto the basis of eigenvectors enables equilibrium states to be preserved [49]. This upwinding is often applied to bathymetry source terms, but we can apply it more generally when derivatives of state variables are involved. We therefore decide to discretise the term $\Delta S_{H_s,i}$ in (3.4), in the latter way:

$$ \Delta S_{H_s,i} = S_{H_s,i+\frac{1}{2}}^{-} (U_{i+\frac{1}{2}}^{L}, U_{i+\frac{1}{2}}^{R}) + S_{H_s,i-\frac{1}{2}}^{+} (U_{i-\frac{1}{2}}^{L}, U_{i-\frac{1}{2}}^{R}) + S_{H_s}^{*} (U_{i+\frac{1}{2}}^{L}, U_{i-\frac{1}{2}}^{R}), \quad (3.8) $$

with

$$ S_{H_s,i \pm \frac{1}{2}}^{\mp} (U_{i \pm \frac{1}{2}}^{L}, U_{i \pm \frac{1}{2}}^{R}) = \frac{1}{2} \left[ \tilde{P} (I \mp \tilde{\Lambda}^{-1} |\tilde{\Lambda}|) \tilde{P}^{-1} \right]_{i \pm \frac{1}{2}} \tilde{S}_{H_s,i \pm \frac{1}{2}}, \quad (3.9) $$

$$ \tilde{S}_{H_s,i \pm \frac{1}{2}} = \left( \frac{A^L + A^R}{2} (H_s^R - H_s^L) \right)_{i \pm \frac{1}{2}}, \quad (3.10) $$

$$ S_{H_s}^* (U_{i+\frac{1}{2}}^L, U_{i-\frac{1}{2}}^R) = \left( \frac{A_{i+\frac{1}{2}}^L + A_{i-\frac{1}{2}}^R}{2} (H_{s,i+\frac{1}{2}}^L - H_{s,i-\frac{1}{2}}^R) \right), \quad (3.11) $$

where $S_{H_s}^*$ is designed to maintain equilibrium states at rest when a high-order reconstruction of state variables is involved. On the other hand, the term $\bar{S}_{\Phi}$ in (3.4) is discretised following [12], by integrating $\Phi$ using the piecewise linear variation underlying the finite element method used in the elliptic step (cf. section §3.2), namely:

$$ \bar{S}_{\Phi} = \frac{1}{8} \left( \begin{array}{c} 0 \\ A_{i+\frac{1}{2}}^{L} (\Phi_{i+1} + 3\Phi_i) + A_{i-\frac{1}{2}}^{R} (3\Phi_i + \Phi_{i-1}) \end{array} \right). \quad (3.12) $$
### 3.1.2. Energy conservative scheme (EC)

We also consider the scheme proposed by Fjordholm et al. (2011), which, in absence of dispersive terms, satisfies a discrete form of the energy equation. Such schemes are suitable for long time computations of smooth solutions. We have adapted the expression of [42] so that the scheme is energy conservative for the Saint-Venant equations expressed in terms of cross-section. The derivation of the numerical entropy conservative flux $F_{i+\frac{1}{2}}$ satisfying Tadmor's condition [50] is written in Appendix B and is expressed as follows:

$$F_{i+\frac{1}{2}} = \begin{pmatrix} \bar{A}_{i+\frac{1}{2}} \bar{u}_{i+\frac{1}{2}} \\ \bar{A}_{i+\frac{1}{2}} \bar{u}_{i+\frac{1}{2}}^2 + g \bar{K}_{i+\frac{1}{2}} \end{pmatrix}, \qquad (3.13)$$

where $\bar{u}_{i+\frac{1}{2}}$, $\bar{K}_{i+\frac{1}{2}}$ are computed using arithmetic mean while we showed that $\bar{A}_{i+\frac{1}{2}}$ has to be an exact mean value:

$$\bar{A}_{i+\frac{1}{2}} \doteq \frac{1}{h_{i+\frac{1}{2}}^R - h_{i+\frac{1}{2}}^L} \int_{h_{i+\frac{1}{2}}^L}^{h_{i+\frac{1}{2}}^R} A(h) \, dh. \qquad (3.14)$$

The source term involving the derivative of the hydrostatic head is now equal to:

$$\Delta S_{H_s,i} = S_{H_s,i+\frac{1}{2}} - S_{H_s,i-\frac{1}{2}}, \qquad (3.15)$$

with

$$S_{H_s,i+\frac{1}{2}} = \begin{pmatrix} 0 \\ \frac{1}{2} \bar{A}_{i+\frac{1}{2}} (H_{s,i+1} - H_{s,i}) \end{pmatrix}. \qquad (3.16)$$

Finally, the $\bar{S}_{\Phi}$ source term is the same as the one used with the Roe flux. This scheme is naturally second-order accurate approximation of (3.1) without using variables reconstructions.
This choice allows to study the impact of the numerical dissipation embedded in the numerical fluxes used in the hyperbolic step. For completeness we stress that for dispersive equations also admitting an exact energy conservation equation, which is not the case of (3.1), the construction of discretely energy conservative schemes is quite more complex. For a particular family of Boussinesq equations, examples are given e.g. in [10, 51].

## 3.2. Elliptic step - Finite Element Method

We are now interested in solving the elliptic equation (2.25) to obtain the auxiliary variable $\Phi$. This is done using a finite element method with $\mathbb{P}_1$ Lagrange elements. With regard to the boundary conditions, we have chosen to impose $\Phi$ such that the Saint-Venant equations are recovered at the edge of the domain $\Omega$. The problem we wish to solve is therefore written as:

$$ (1 + \mathcal{T}) \Phi = -H_x, \quad \text{in } \Omega, \quad (3.17) $$

$$ \Phi = -H_x, \quad \text{on } \partial\Omega, \quad (3.18) $$

with $\mathcal{T}$ and $H$ defined by (2.17) and (2.23), respectively.

The variational formulation of this equation is obtained by multiplying by a test function $v$ and integrating over the domain $\Omega$:

$$ \int_{\Omega} (1 + \alpha) \Phi v \, dx + \int_{\Omega} \beta^* \Phi_x v \, dx - \int_{\Omega} (\gamma^* \Phi_x)_x v \, dx = - \int_{\Omega} H_x v \, dx. \quad (3.19) $$

Differently from (3.3) in which $U_i$ denote cell averages, our main unknown is now the array $\{\Phi_i\}_{1\le i\le N}$ which contains approximations of the point-wise values of the dispersive source $\Phi$ at the mesh nodes $(x_i)_{1\le i\le N}$. We use a classical piecewise linear $\mathbb{P}_1$ finite elements expansion $v_h$ defined by as

$$ v_h(x) = \sum_{j=1}^{N} v_j \varphi_j, \quad (3.20) $$
where the basis functions $\varphi_j(x)$ are the usual piecewise linear Lagrange polynomials. It is important to note that the $i$th node of the FE discretisation coincide with the cell centre of cell $i$ in the FV discretisation. Substituting (3.20) in (3.19) leads to the following system of linear equations:

$$ \sum_{j=1}^{N} a(\varphi_j, \varphi_i) \Phi_j = b(\varphi_i), \qquad (3.21) $$

with

$$ \begin{aligned} a(\varphi_j, \varphi_i) &\doteq \int_{\Omega} (1+\alpha) \varphi_j \varphi_i \,dx + \int_{\Omega} \beta^* \varphi_{j,x} \varphi_i \,dx + \int_{\Omega} \gamma^* \varphi_{j,x} \varphi_{i,x} \,dx, \\ b(\varphi_i) &\doteq -\int_{\Omega} H_x \varphi_i \,dx + (\gamma^* \Phi_x \varphi_i)|_{\partial \Omega}. \end{aligned} $$

Calculating the matrix of components $a(\varphi_j, \varphi_i)$ and the right-hand side $b(\varphi_i)$, only require the estimation of the volumic integrals. Here we use a Gauss second-order quadrature method. It is important to note that the mass matrix only depends on $\alpha$, $\beta^*$ and $\gamma^*$. Since these coefficients are independent on time (we assume the channel to remain steady), the coefficients of the linear system are estimated once for ever on the first time step. The resulting tridiagonal system of equations is solved using the Thomas algorithm [52].

### 3.3. Time scheme and friction

#### 3.3.1. Time integration

For time integration of (3.4), we use the method of Heun [53] and SSP Runge-Kutta 3 (RK3) [54] to obtain second and third-order schemes, respectively. Both are predictor-corrector methods, where Heun's method consists
of 2 steps which can be written as follows:

$$ 
\begin{aligned}
\boldsymbol{U}^p &= \boldsymbol{U}^n + \Delta t \boldsymbol{\Psi} (\boldsymbol{U}^n), \\
\boldsymbol{U}^{2p} &= \boldsymbol{U}^p + \Delta t \boldsymbol{\Psi} (\boldsymbol{U}^p), \\
\boldsymbol{U}^{n+1} &= \frac{1}{2} (\boldsymbol{U}^n + \boldsymbol{U}^{2p}),
\end{aligned} 
\qquad (3.22) $$

while RK3 involves 3 steps:

$$ 
\begin{aligned}
\boldsymbol{U}^p &= \boldsymbol{U}^n + \Delta t \boldsymbol{\Psi} (\boldsymbol{U}^n), \\
\boldsymbol{U}^{2p} &= \frac{3}{4} \boldsymbol{U}^n + \frac{1}{4} (\boldsymbol{U}^p + \Delta t \boldsymbol{\Psi} (\boldsymbol{U}^p)), \\
\boldsymbol{U}^{n+1} &= \frac{1}{3} \boldsymbol{U}^n + \frac{2}{3} (\boldsymbol{U}^{2p} + \Delta t \boldsymbol{\Psi} (\boldsymbol{U}^{2p})),
\end{aligned} 
\qquad (3.23) $$

with $\Delta t^{n+1} = t^{n+1} - t^n$ and $\boldsymbol{\Psi}$ the RHS of (3.4):

$$ \boldsymbol{\Psi} = -\frac{1}{\Delta x} (\boldsymbol{F}_{i+\frac{1}{2}} - \boldsymbol{F}_{i-\frac{1}{2}}) + \frac{1}{\Delta x} \Delta \boldsymbol{S}_{H_s, i} + \bar{\boldsymbol{S}}_{\Phi}. \qquad (3.24) $$

On the other hand, $\Delta t$ is computed at each time step and is based on the Courant–Friedrich–Lewy (CFL) condition:

$$ \Delta t^{n+1} = \text{CFL} \frac{\Delta x}{\max (|u_i^n| + \sqrt{g\bar{h}_i^n})}. \qquad (3.25) $$

In all numerical computations we use a CFL number of 0.25, unless otherwise specified.

### 3.3.2. Friction terms: semi-implicit integration

To take account of bottom friction, an additional source term is explicitly added to (3.1), which is modified as

$$ \boldsymbol{S} = \begin{pmatrix} 0 \\ A\Phi + AH_{s,x} - gAS_f \end{pmatrix}, $$
where $S_f$ is evaluated using the Gauckler–Mannig–Strickler empirical formula:

$$S_f = \frac{(Au)^2}{K_s^2 A^2 R_h^{4/3}}, \qquad (3.26)$$

with $K_s$ the Strickler coefficient and $R_h$ the hydraulic radius. This friction term is treated numerically in a semi-implicit way, avoiding numerical oscillations as reported in [55]. This scheme consists of first estimating the $U^*$ solution without the friction term:

$$U_i^* = U_i^n - \Delta t \Psi.$$

Then friction is added as follows:

$$A_i^{n+1} = A_i^*, \qquad (3.27)$$

$$(Au)_i^{n+1} = \frac{(Au)_i^*}{1 + \Delta t g \frac{|(Au)_i^n|}{K_s^2 A_i^* (R_h^{4/3})_i^n}} \qquad (3.28)$$

The FE and FV parts of the solver were verified independently using the Method of Manufactured Solution (MMS) [56],[57]. This method is based on imposing an exact source term calculated from an analytical solution. In this case, for the hyperbolic part we considered the Serre soliton [58], then deduced the right-hand side member which we then imposed as an exact source term. On the other hand, and in the same way for the elliptic part, we checked with two forms of solution for $\Phi$, namely a sine and a gaussian function with non constant values for $\alpha, \beta$ and $\gamma$. In both cases we were able to find the order of accuracy of the schemes involved.
## 4. Analytical solution in prismatic channel of trapezoidal cross-section

### 4.1. Derivation

Winckler and Liu derived an analytical solution to their model for prismatic channels of arbitrary cross-section. However, in addition to approximating the continuity equation, they also made other approximations in their derivation, following Mei's [59] derivation method. With our modified system we can derive a solitary wave solution for the special case of prismatic channels with symmetrical trapezoidal cross-sections. The non-symmetric case can be derived using the same method and is not treated in this work. We consider the system (2.30) and seek for a solution that propagates at celerity $c$, thus introduce a moving coordinate $\xi = x - ct$. For a solitary wave we prescribe $u = 0$ and $\eta = 0$ for $\xi = \pm\infty$. The continuity equation can then be integrated as:

$$u = \frac{c}{\epsilon} \left( 1 - \frac{A_0}{A} \right). \qquad (4.1)$$

We proceed in the same way by applying this change of coordinates to the momentum equation. Using (4.1) to remove $u$, and integrating once leads to:

$$\frac{1}{2} \left( \frac{A_0^2}{A^2} - 1 \right) + \frac{\epsilon g}{c^2} \eta - \mu^2 \gamma^* \left( \frac{A_0}{A} \right)'' = 0, \qquad (4.2)$$

where the prime is the derivative with respect to $\xi$. We now multiply by $(\frac{A}{A_0})'$ and integrate once more. The following equation is obtained after some manipulations:

$$3\mu^2\gamma^* \left(\frac{A_0}{A}\right)'^2 = \left(\frac{A_0}{A}\right)^3 - 3\left(\frac{A_0}{A}\right) + 2 + 6\frac{\epsilon g}{c^2} \left(\frac{A_0}{A}\eta - \int \frac{A_0}{A} d\eta\right). \quad (4.3)$$
<figure><img src="image_3.png" /><figcaption>Figure 2: Notation used for symmetric trapezoidal cross-section</figcaption></figure>

Noting that $A' = B\eta'$ with B the width at $z = \eta$, and multiplying by $(\frac{A}{A_0})^4$ we obtain:

$$3\mu^2\gamma^*\frac{B^2}{A_0^2}\eta'^2 = \left(2\frac{A}{A_0}+1\right)\left(\frac{A}{A_0}-1\right)^2\frac{A}{A_0}-6\frac{\epsilon g}{c^2}\left(\frac{A}{A_0}\right)^4\left(\int\frac{A_0}{A}d\eta-\frac{A_0}{A}\eta\right). \quad (4.4)$$

For a symmetrical trapezoidal cross-section with bank slope tan $\theta$, bottom width b and width-averaged water depth above the bottom h (Figure 2), we introduce $\tau(h) \doteq \frac{h}{b} \cot \theta$ that represents the shape of the wetted cross-section A. We thus have:

$$A = \left[ 1 + \tau_0 + (1 + 2\tau_0) \frac{\eta}{h_0} + \tau_0 \frac{\eta^2}{h_0^2} \right] bh_0, \quad (4.5)$$

$$A_0 = (1 + \tau_0) bh_0, \quad (4.6)$$

$$\frac{A}{A_0} = 1 + \frac{1 + 2\tau_0}{1 + \tau_0} \frac{\eta}{h_0} + \frac{\tau_0}{1 + \tau_0} \frac{\eta^2}{h_0^2} = \frac{1}{1 + \tau_0} \frac{h}{h_0} \left( 1 + \tau_0 \frac{h}{h_0} \right), \quad (4.7)$$

$$B = A_{\eta} = \left( 1 + 2\tau_{0}\frac{h}{h_{0}} \right) b, \quad (4.8)$$

where $\tau_0 \doteq \tau(h_0)$. Substituting in (4.4) and performing simple algebraic manipulations leads to:

$$3\mu^2 (1 + \tau_0) \frac{\gamma^*}{h_0^2} \left(1 + 2\tau_0 \frac{h}{h_0}\right)^2 \eta'^2 = \varphi\left(\frac{h}{h_0}\right) - \frac{6\epsilon}{F_0^2} \psi\left(\frac{h}{h_0}\right), \quad (4.9)$$
where $F_0 \doteq \frac{c}{\sqrt{gh_0}}$ is a Froude number, and:

$$ \varphi(x) \doteq \left( \frac{2}{1 + \tau_0} x (1 + \tau_0 x) + 1 \right) (x (1 + \tau_0 x) - 1 - \tau_0)^2 x (1 + \tau_0 x), \quad (4.10) $$

$$ \psi(x) \doteq x^4 (1 + \tau_0 x)^3 \left( (1 + \tau_0 x) \ln \frac{x (1 + \tau_0)}{1 + \tau_0 x} + \frac{1}{x} - 1 \right). \quad (4.11) $$

At the solitary wave crest, i.e. $\frac{h}{h_0} = \frac{h_{max}}{h_0} = \frac{h_0+a_0}{h_0} = 1 + \epsilon$, the slope of the free surface elevation is zero, $\eta' = 0$. We can then express the celerity of the solitary wave $c$ as a function of the maximum water height $1 + \epsilon$:

$$ F_0^2 = \frac{c^2}{gh_0} = 6\epsilon \frac{\psi(1 + \epsilon)}{\varphi(1 + \epsilon)}. \quad (4.12) $$

Substituting this in (4.9), we finally obtain the following closed form solution:

$$ \frac{x - ct}{h_0} = \mu \sqrt{3(1 + \tau_0) \frac{\gamma^*}{h_0^2}} \int_{h/h_0}^{1+\epsilon} \frac{1 + 2\tau_0 H}{\sqrt{\varphi(H) - \frac{\varphi(1 + \epsilon)}{\psi(1 + \epsilon)} \psi(H)}} \, dH, \quad (4.13) $$

where $H$ is a dummy variable. This integral can be computed numerically using quadrature methods. In this case, we consider this solution for $x - ct < 0$, because time is passing. For a prismatic channel of rectangular cross-section $\tau_0 = 0$ and $\gamma_{rect}^* \doteq \frac{1}{3}h_0^2$ [39], the solution (4.13) leads to:

$$ \frac{x - ct}{h_0} = \mu \int_{h/h_0}^{1+\epsilon} \frac{dH}{\sqrt{\varphi(H) - \frac{\varphi(1 + \epsilon)}{\psi(1 + \epsilon)} \psi(H)}}. \quad (4.14) $$

Note that the integral in (4.13) has a singularity at $\frac{h}{h_0} = 1 + \epsilon$, where the function under the square root of the integral is zero. We show in Appendix C that this integral converges in the neighbourhood of this singularity.
## 4.2. Comparison of rectangular and trapezoidal cross-sections

The solution (4.13) is plotted in Figure (3a) for $\tau_0 = 0, \frac{1}{2}, 1, \frac{3}{2}$, with the same shape factor $\frac{b}{h_0} = 2$. Except for the rectangular case, where $\gamma^*$ has an analytical expression, we solve the BVP (2.13) using the PDE solver FreeFem++, for $\tau_0 \neq 0$. We see that the effective wavelength depends on $\tau_0$. For small values of $\tau_0$, e.g. $\tau_0 = \frac{1}{2}$, it is shorter than in the rectangular case despite $A_0$ being larger. To investigate this behavior, $\gamma^*$ is plotted as a function of $\tau_0$ for different aspect ratios $\frac{b}{h_0}$. We see that $\gamma^*$ is smaller than $\gamma_{rect}^*$ for $\tau_0 < 0.57$ and $\frac{b}{h_0} = 2$. The greater $\frac{b}{h_0}$, the faster $\gamma^*$ moves away from $\gamma_{rect}^*$. On the other hand, for the same geometric configurations, we look at the dimensionless celerity of the solitary wave, i.e. the Froude number $F_0$, as a function of its dimensionless maximum amplitude $\epsilon$, given by (4.12) as in figure 3b. Firstly, we can see that the celerity $c$ is greater in a rectangular channel than in a trapezoidal channel, for any $\epsilon$. In other words, the larger the cross-section is, the slower the wave, as already shown in [38].

## 4.3. Numerical verification of the coupling

The solitary wave solution (4.13) is now used to verify the whole scheme. This test case is relevant for studying the stability and conservation properties of our numerical scheme. We study the soliton propagation in a prismatic channel ($\alpha = \beta = 0$ and $\gamma = \text{cst}$) of trapezoidal cross-section with $h_0 = b = 1m$ and $\tan\theta = \frac{1}{3}$, i.e. $\tau_0 = 3$. For this geometry, we solve the BVP (2.13) using FreeFem++, and obtain $\gamma^* = 0.6406m^2$. The solitary wave is initialized with $\epsilon = 0.2$ in a domain $x \in [-10.5\lambda, 10.5\lambda]$, where $\lambda = 19m$ is defined as the distance to the wave peak where the free surface elevation is 1% of the wave amplitude. As boundary conditions we prescribe
<figure><img src="image_4.png" /><figcaption>Figure 3: (a) Trapezoidal channel cross-section: Analytical solution 4.13 for various values of $\tau_0$, (b): soliton speed as a function of $\epsilon$, for the same values of $\tau_0$. (c): ratio of $\gamma^*/\gamma_{rect}^*$ as a function of $\tau_0$, for various values of the shape factor $\frac{b}{h_0}$.</figcaption></figure>
the analytical solution, i.e. $u(-10.5\lambda, t) = u_{exact}$ and $h(x = 10.5\lambda, t) = h_{exact}$.
A convergence study is shown in figure 4a when the soliton has travelled
$4.24\lambda$; for $\Delta x = 0.8, 0.4, 0.2, 0.1, 0.05, 0.025, 0.0125$ m corresponding to a
number of cells per wavelength $N_\lambda \approx 24, 48, 95, 190, 380, 760, 1520$. We
look at the free surface elevation L2-norm error $\| e \|_{L^2}$ as a function of the
number of cells $N = \frac{L}{\Delta x}$, where $L$ is the length of the domain. We plotted
the convergence of a first (Roe O1), second (Roe O2) and third (Roe O3)
order Roe scheme, in order to show the benefits of a higher order. For each
scheme we see that we are close to the theoretical slopes. Furthermore, we
can see that the errors given by the second-order scheme with the Roe flux
give slightly smaller errors than with EC for the same grid size and Roe O3
gives a convergence speed of order 2 because of the finite element solver.
However, this resolution reduces the error by a factor of 2 compared with
Roe O2. As expected, the errors of the first-order scheme are greater than
those of higher-order schemes for all grid sizes. Figure 4b shows the shapes
of the solutions at $t = \frac{4.24\lambda}{c}$ for $\Delta x = 0.1$, and it can be seen that the first-
order scheme dissipates a lot, impacting both the shape and the amplitude
of the soliton. Conversely, the results given for the higher-order schemes give
a solution very close to the analytical solution. Following this verification of
the schemes, we will only consider the second and third-order Roe and EC
schemes, as the numerical dissipation with a first-order scheme will be too
high for long propagation times and reasonable mesh sizes.
<figure><img src="image_5.png" /><figcaption>Figure 4: (a) Solitary wave case. Grid convergence of the solitary wave case for all schemes; (b) Shape of the free surface for= all schemes at $t = \frac{4.24\lambda}{c}$. The wave propagates from left to right.</figcaption></figure>

## 5. Application to the study of Favre waves

In this section we will now apply the numerical schemes to the case of Favre wave propagation. To this end, we represent a rapid gate opening on water at rest ($u_0 = 0 m.s^{-1}$), by initializing a jump from the Rankine-Hugoniot relations:

$$u_0 - c_b = \sqrt{g \frac{A_1 K_1 - K_0}{A_0 (A_1 - A_0)}}, \quad (5.1)$$

$$u_1 - c_b = \sqrt{g \frac{A_0 K_1 - K_0}{A_1 (A_1 - A_0)}}, \quad (5.2)$$

where $A(h) = bh(1 + \tau(h))$ and $K(h) = bh^2(\frac{1}{2} + \frac{1}{3}\tau(h))$, for a trapezoidal section; and where indices 0 and 1 represent the variables upstream and downstream of the jump, respectively, as shown in figure 5. In the following, initial conditions are smoothed from the relationships (5.1), (5.2):
<figure><img src="image_6.png" /><figcaption>Figure 5: Notation of variables of interest for Favre waves propagation</figcaption></figure>

$$h(x, t = 0) = \frac{h_1 - h_0}{2} \tanh \left( \frac{x - x_c}{a} \right) + h_1 - \frac{1}{2}\epsilon h_0, \quad (5.3)$$

$$u(x, t = 0) = \frac{u_1 - u_0}{2} \tanh \left( \frac{x - x_c}{a} \right) + \frac{1}{2} u_1, \quad (5.4)$$

where $x_c$ is the abscissa where the jump is centered and $a = 0.5$ m, the stiffness of the initial jump. For the boundary conditions, we set $u(x = 0, t) = 0$ and $h(x = L, t) = h_1$. We look at how the maximum amplitude $a_{max}$ of the first wave and the wavelength $\lambda$ measured between the first two waves evolve for long propagation times. We then perform a validation from the experiments performed by Favre [28] and Treske [30].

### 5.1. Long propagation time behaviour

To investigate the influence of the numerical scheme on the behavior of Favre waves for long propagation times, we initialize a jump at $x_c = 900$ m in a domain $x \in [0, 1000]$ m, the wave propagating towards $x = 0$. We consider a rectangular channel ($\tau_0 = 0$) with $h_0 = 0.2$, for three initial jump amplitudes $h_1$ corresponding to $\epsilon = \frac{h_1 - h_0}{h_0} = 0.1, 0.2, 0.3$, i.e. for a
Froude number computed in a frame of reference moving with the jump
Fr $\doteq \frac{|u_0-c_b|}{\sqrt{gh_0}} \approx 1.07, 1.15, 1.22$, respectively. For spatial discretization, we
considered 3 meshes, with $\Delta x = 0.1, 0.05, 0.025 \, m$. Figure 6 shows the evo-
lution of $\lambda$ and $a_{max}$ as a function of a dimensionless time $t' = \sqrt{\frac{g}{h_0}}t$, for each
configuration. In addition to the numerical results obtained with the second
and third-order schemes of Roe and EC, we plot the logarithmic relation of $\lambda$
given in [27] (eq. (5.28)) and established from the Korteweg-de Vries equa-
tion. Firstly, we see that for $\epsilon = 0.1$, the EC scheme converges more rapidly
with mesh refinement, in particular that the difference between $\Delta x = 0.05 \, m$
and $0.025 \, m$ is smaller than that obtained with Roe's solver. In addition, EC
seems to slightly overestimate $\lambda$ for $\Delta x = 0.1, 0.05 \, m$, whereas the opposite
is observed for Roe's solver. Furthermore, we observe that the larger $\epsilon$ is, the
more difficult it is to converge in mesh, where in particular calculations blow
up with EC for $\Delta x = 0.025 \, m$ and $\epsilon = 0.2, 0.3$, which explains why these
cases are not depicted on the figures. This can be explained by the fact that
the velocities are higher, increasing instabilities due to the centered nature
of EC. Note also that time-explicit schemes such as Heun and RK3 produce
energy that can lead to instabilities, since there is no artificial dissipation
introduced with EC. Also, for EC and with the thinnest meshes, we can ob-
serve a logarithmic evolution of the wavelength as a function of time, insofar
as we obtain straight lines with different slopes to those given by the [27]
relationship, but this can be explained by the fact that we are not working
with the same equations.

In the absence of friction, the wavelength increases continuously with time,
while the maximum amplitude seems to tend towards a horizontal asymptote
for any $\epsilon$, as can be seen in the graphs on the right of the same figure. Similarly to the wavelength, we can see that numerical dissipation is significant with the Roe's scheme for all grid sizes, while with EC and $\Delta x = 0.1 \, m$ we obtain a correct estimate, particularly for long times, in comparison with solutions obtained with thinner meshes. Although 3rd-order for the hyperbolic part (Roe O3) reduces the numerical dissipation, we can see that it remains significant, especially as $\epsilon$ increases. For long times, upwind schemes of higher order of accuracy (4th, 5th) should be used to obtain results similar to EC.

The EC scheme seems to better represent the propagation physics of the equations considered in the sense that, in the absence of friction, no artificial dissipation is introduced. For realistic applications at large scales, however, bottom friction dissipation can be significant, depending of course on the materials constituting the channel. Therefore we will now look at, how bottom friction, influences the propagation dynamics represented by the numerical schemes. We consider the same initial configuration as above and a Strickler coefficient $K_s = 100m^{1/3}s^{-1}$. Figure 7 shows once again the evolution of $\lambda$ and $a_{max}$ as a function of $t'$ for each $\epsilon$. Taking friction into account gives a different dynamic for the wavelength with the Roe's schemes. For long times, it stabilises and then increases again, as for example for Roe O2 with $\Delta x = 0.05 \, m$ and $\epsilon = 0.2$. Concerning the evolution of $\lambda$ with the EC scheme, we can see that the wavelength seems to grow logarithmically until about $t' = 1000$. Thereafter, it seems to stabilize, also tending towards a horizontal asymptote. On the other hand, for the maximum amplitude of the first wave, the friction term rapidly affects its evolution, i.e. the larger $\epsilon$ is,
the sooner $a_{max}$ reaches its maximum before decreasing in an almost linear way. We can see that the estimation of this amplitude is particularly accurate with $\Delta x = 0.1$ m, even for short times ($t' < 500$), compared to the case without friction. Compared to the frictionless case, Roe O3 and EC give similar results, particularly for the maximum amplitude of the first wave. However, for the wavelength, we can see that Roe O3 still requires greater spatial discretisation.

These results confirm our claims concerning the choice of decentered numerical fluxes for modelling such waves at long propagation times. Indeed, as we can see from the results given with Roe's scheme, the evolution of the quantities of interest are poorly represented, largely dominated by numerical dissipation. Conversely, the EC flux gives satisfactory results, although clearly some additional care is required in the time stepping to ensure their stability.

## 5.2. Comparison with experiments of Favre (1935) and Treske (1994) in a rectangular channel

We now apply Winckler and Liu's reformulated model to the case of Favre's channel [28], which is 73.58 meters long, 0.42 meters wide and 0.4 meters high. We still consider the case of a rapid opening of the gate, with the initialization (5.3) where the jump is now centered at $x_c = 68.58$ m and the water height before opening is $h_0 = 0.2$ m. Favre carried out various tests, varying the flow rate at $x = 73.58$ m, resulting in a mobile hydraulic jump with different water heights $h_1$, i.e. different Froude numbers Fr. For each Froude number tested, he measured the wavelength $\lambda$ and the maximum amplitude of the first wave $a_{max}$ after approximately 63.5 m of trav-
<figure><img src="image_7.png" /><figcaption>Figure 6: Favre waves case in a rectangular cross-section channel. $\lambda$ and $a_{max}$, as functions of $t'$, given for all schemes and different grid sizes (no friction). Wavelength relation from [27] is given. Top: $\epsilon = 0.1$, Middle: $\epsilon = 0.2$, Bottom: $\epsilon = 0.3$</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Favre waves case in a rectangular cross-section channel. $\lambda$ and $a_{max}$, as functions of $t'$, given for all schemes and different grid sizes (with friction). Wavelength relation from [27] is given. Top: $\epsilon = 0.1$, Middle: $\epsilon = 0.2$, Bottom: $\epsilon = 0.3$.</figcaption></figure>
elling distance. We tested various jump amplitudes ranging from $\epsilon = 0.02$ to $\epsilon = 0.42$ with a step of 0.04. For data comparison, we stopped the simulation when the maximum of the first wave arrived at $x = 10m$. As in the previous section, we look at the results given with Roe and EC fluxes, but with and without the friction term, where $K_s = 100m^{1/3}s^{-1}$. The figure (8) shows the evolution of $\lambda$ and $a_{max}$ as a function of Fr. Our numerical results reached convergence with $\Delta x = 0.009m$. In addition, we have also added data from Treske's experiments in a rectangular channel. For each case, we observed a correct decay of the wavelength as the Froude increases. However, a larger overestimation of the wavelength is obtained with increasing Fr, which has also been observed with the SGN equations in [43]. This overestimation is more pronounced when compared with Treske's data; note that both Favre and Treske, did not consider the increase of the wavelength with time. Treske does not explicitly state at which position he measured $\lambda$ and $a_{max}$ in his 100-meter-long channel. It is therefore difficult to compare the results. Furthermore, and in line with our previous remarks, we can see that the friction term had almost no influence on the wavelength for all schemes, which can be explained by the short propagation time of the wave in this channel. Conversely, we can see that we have good estimates of maximum amplitudes, with all schemes. In particular, the effect of friction is more noticeable here, where for $Fr < 1.25$, the estimate is very satisfactory. To be able to estimate more precisely for $Fr > 1.25$, a wave breaking model would have to be introduced, since as specified in [30] for $Fr > 1.3$ the first wave starts to break. To summarize, the numerical model, both with the Roe's and EC scheme, gives promising results, although comparison with [28] and
<figure><img src="image_9.png" /><figcaption>Figure 8: Favre waves in a rectangular cross-section channel. Comparison of (a) $\lambda$ and (b) $a_{max}$ as functions of Fr, for all schemes, with and without friction; and with Favre (1935) and Treske (1994)'s data performed in a rectangular channel.</figcaption></figure>

[30] remains cautious.

### 5.3. Comparison with experiments of Treske (1994) in a trapezoidal channel

We now extend our validation of the model (2.30) by comparing with Treske's tests in a trapezoidal cross-section channel. Treske [30] carried out measurements in a channel whose dimensions are: $L = 130 \text{ m}$, $b = 1.24 \text{ m}$ and $\tan \theta = \frac{1}{3}$, with an initial water height $h_0 = 0.16 \text{ m}$, i.e. $\tau_0 \approx 0.387$, for which we obtained $\gamma^* = 0.02477 \text{ m}^2$. With these tests, Treske highlighted a transition when $Fr > 1.15$, the propagation dynamics deviating from those observed in a rectangular channel, with refractions and breaking on the banks giving way to a fishtail pattern downstream of the wave front. Following this transition, and before the second one (the fully breaking bore), he was able to measure a drop in wavelength at the center of the channel. In addition to the measurements at the center, he also measured on the banks, where he
observed greater amplitudes on the banks than at the center of the channel. As mentioned in the previous section, it is not clear at which position he measures $a_{max}$ and $\lambda$, as he presents measurements of the free surface as a function of time at two different positions $x = 72m$ and $x = 100m$. We choose to initialize our simulations with (5.3) and $x_c = 125m$, then stop our simulations when the maximum of the first wave is at $x = 72m$. As in the previous section, we compared the numerical results given with Roe's and EC schemes, with and without friction. For each case, we discretized in space with $\Delta x = 0.01625m$. In addition to the numerical results, we plotted on figure 9 the spanwise-weighted average of Treske's data, assuming that the profile between banks and channel center is parabolic. With regard to wavelength, our simulations once again captured the decrease in wavelength as $Fr$ increases, notably for $Fr < 1.15$, before the transition. Beyond this threshold, wavelength estimation is rather poor, due in particular to the short experimental wavelength at the center of the channel, a phenomenon that cannot be reproduced with our 1D model. Conversely, for $a_{max}$, we see that the estimation is much better, especially for $Fr < 1.15$, even slightly beyond which the deviation remains acceptable. Moreover, we observe that friction still has no significant role on $\lambda$, whereas it is even more significant than in the rectangular case, for $a_{max}$ with great $\epsilon$. This may be explained by the fact that the hydraulic radius $R_h$ is smaller than in Favre experiments. Although it is once again difficult to estimate the wavelength adequately, we can see that despite the 1D modelling the model manages to reproduce the main features for the maximum amplitude of the first wave.
<figure><img src="image_10.png" /><figcaption>Figure 9: Favre waves in a trapezoidal cross-section channel. Comparison of (a) $\lambda$ and (b) $a_{max}$ as functions of Fr for all schemes, with and without friction; and with Treske (1994)'s data performed in a trapezoidal channel.</figcaption></figure>

## 6. Conclusion

This study proposes a reformulation of the Winckler and Liu model, a Boussinesq-type model adapted to wave propagation in channels with variable cross-sections. In particular, this new model makes adding the dispersive term to an existing code that solves the Saint-Venant equations easy. The hyperbolic and dispersive parts of the problem are approached using a well-balanced finite volume and finite element method. A solitary wave solution without additional approximations is derived, allowing the numerical verification to be carried out accurately. Different numerical flux formulations for the hyperbolic part are studied through the case of Favre waves propagations, namely the traditional and often used Roe formulation and an energy conservative flux adapted to the section-averaged shallow water equations. Thus, through long propagation times and in the absence of friction, we show
the need to use high orders of accuracy for flux such as Roe in comparison
with EC, to attenuate numerical dissipation. The addition of friction shows
the regularisation that such artificial dissipation can have, which must there-
fore be considered in the approximation of Boussinesq-type equations. On
the other hand, by validating the proposed model with data from Favre and
Treske's laboratory experiments, we show that the model can represent the
key phenomena of propagation, especially for the amplitudes, even when the
dynamics is three-dimensional. The approximations made on the initialisa-
tion of the test cases and the lack of details on the experimental set-up, lead
to differences in the wavelength, We highlight the importance of taking into
account the increase in this parameter to arrive at accurate comparisons.
This promising reformulation will be the subject of future validation with
new laboratory tests on Favre waves focusing on the time evolution of the
characteristic quantities, maximum amplitude and characteristic wavelength,
studied in the present work.

7. Acknowledgements

Authors express their gratitude for the financial support provided by the
French ANRT (Association Nationale de la Recherche et de la Technologie)
with the CIFRE Grant 2021/1094 towards Bastien Jouy's Ph.D. research
program.
## Appendix A. Solutions of $\alpha$, $\beta$ et $\gamma$ for triangular cross-section

We show in section 2.4 that the coefficient $\gamma^*$ is positive (i.e $\gamma$ is negative), regardless the geometry of the cross-section. However, Winckler and Liu present in their equation (3.57) in [39] a coefficient $\gamma$ that can be either positive or negative depending on the width of the unperturbed free surface $B_0$. By doing their derivation for $\alpha$, $\beta$ and $\gamma$ for a symmetric triangular cross-section we found different expressions. We found a mistake in their average over the cross-section where the bottom is now depending of $y$ (compared to rectangular cross-section):

$$z = -h_0(x, y) = \frac{2H_0(x)}{B_0(x)}|y| - H_0(x),$$

where $H_0(x)$ is the maximum depth at rest of the triangular cross-section. For this geometry the solution to the BVP (2.13) is as expressed in [39]:

$$\chi_2(y, z) = -\frac{1}{4}(2H_0z + z^2 + y^2).$$

Then as an example, we show next our derivation of $\gamma$ (2.9), for a triangular cross-section, first we need $\tilde{\chi}_2$:

$$\tilde{\chi}_2 = \frac{1}{B_0} \int_{-B_0/2}^{B_0/2} \chi_2 \, dy = -\frac{B_0^2}{48}, \qquad (A.1)$$
and for $\overline{\chi_2}$, using the following change of variable $y = \mp \frac{B_0(h_0-H_0)}{2H_0}$ and $A_0 = \frac{1}{2}B_0H_0$, we obtain:

$$ 
\begin{aligned}
\overline{\chi_2} &= \frac{1}{A_0} \int_{-B_0/2-h_0}^{B_0/2} \int_{-h_0}^{0} \chi_2 \, dy \, dz, \\
&= \frac{2}{H_0^2} \int_{0}^{H_0} \left( \int_{-h_0}^{0} \chi_2 \, dz \right) \, dh_0, \\
&= \frac{H_0^2}{8} - \frac{B_0^2}{48}. \qquad (A.2)
\end{aligned}
$$

Finally from (A.1) and (A.2), the expression of $\gamma^*$ is:

$$ \gamma^* = \frac{H_0^2}{8} + \frac{B_0^2}{96}. \qquad (A.3) $$

Averages over the cross-section are also needed both for $\alpha$ (2.7) and $\beta$ (2.8). The coefficients computed in the same way as $\gamma^*$ are expressed as follows, for a triangular cross-section:

$$ 
\begin{aligned}
\alpha &= \left( \frac{B_0}{96} - \frac{H_0^2}{8B_0} \right) B_{0,xx} + \left( \frac{1}{32} + \frac{H_0^2}{8B_0^2} \right) (B_{0,x})^2 + \left( -\frac{B_0^2}{96H_0} - \frac{5H_0}{24} \right) H_{0,xx} \\
&\quad + \left( \frac{B_0^2}{96H_0^2} + \frac{1}{24} \right) (H_{0,x})^2 - \left( \frac{B_0}{24H_0} + \frac{H_0}{6B_0} \right) H_{0,x} B_{0,x}, \qquad (A.4)
\end{aligned}
$$

$$ \beta = \left( -\frac{B_0}{32} - \frac{H_0^2}{8B_0} \right) B_{0,x} + \left( -\frac{3H_0}{8} - \frac{B_0^2}{96H_0} \right) H_{0,x}. \qquad (A.5) $$

## Appendix B. Energy Conservative flux for 1D Saint-Venant equations on variable sections

The Energy Conservative flux of Fjordholm et al [42], is expressed for the Saint-Venant equations averaged over the vertical. Following Tadmor's
condition [50], we built an EC scheme adapted to the cross-section averaged Saint-Venant equations. Starting with the latter:

$$A_t + (Au)_x = 0, \qquad (B.1)$$

$$(Au)_t + (Au^2 + gK)_x = 0, \qquad (B.2)$$

where the expression of K is given in (2.28). First, we introduce the kinetic energy $\kappa \doteq \frac{1}{2}u^2$ and we seek an entropy equation, multiplying (B.2) by $u$ and (B.1) by $-\kappa$, we obtain:

$$(A\kappa)_t + (Au\kappa)_x + gAuh_x = 0. \qquad (B.3)$$

We then multiply (B.1) by $gh$ and sum to (B.3), we obtain the following entropy equation using the relation $\partial(Ah) = h\partial A + \partial K$:

$$E_t + F_{E,x} = 0, \qquad (B.4)$$

where $E(U)$ and $F_E(U)$ are respectively an entropy function and an entropy flux defined as follows:

$$E \doteq g(Ah - K) + A\kappa, \qquad (B.5)$$

$$F_E \doteq Au\kappa + ghAu. \qquad (B.6)$$

We also introduce the entropy variable $W(U) = \begin{pmatrix} gh - \kappa \\ u \end{pmatrix}$ which satisfies the following relationships:

$$\partial E = W^t \partial U, \quad \partial F_E = W^t \partial F,$$

with $U$ and $F$ given by (3.2). We can now compute the entropy potential $v$ as defined in [50]:

$$v = W^t F - F_E = guK. \qquad (B.7)$$
Finally, to check the convexity of E we postulate the existence of a function B(h) such that

$$ \partial A = B(h) \partial h, $$

and we compute the following jacobian matrix

$$ M \doteq \frac{\partial U}{\partial W} = \frac{B}{g} \begin{pmatrix} 1 & u \\ u & u^2 + A \frac{g}{B} \end{pmatrix}. \qquad (B.8) $$

Then, for any $V \in \mathbb{R}^2$ we have

$$ VMV = \frac{\sigma}{g} (V_1 + uV_2)^2 + AV_2^2, $$

which is guaranteed positive as soon as the width B is positive.

We now applied to (3.4), without source terms:

$$ \Delta x \frac{d}{dt} U_i + F_{i+\frac{1}{2}} - F_{i-\frac{1}{2}} = 0. \qquad (B.9) $$

We consider classical characterization of entropy conservative fluxes due to [50] which states that a numerical flux is entropy conservative if:

$$ \llbracket v \rrbracket_{i+\frac{1}{2}} = \llbracket W \rrbracket_{i+\frac{1}{2}} F_{i+\frac{1}{2}}, \qquad (B.10) $$

with an EC flux form taken as follows:

$$ F_{i+\frac{1}{2}} = \begin{pmatrix} \bar{A}_{i+\frac{1}{2}} \bar{u}_{i+\frac{1}{2}} \\ \bar{A}_{i+\frac{1}{2}} \bar{u}_{i+\frac{1}{2}}^2 + g \bar{K}_{i+\frac{1}{2}} \end{pmatrix}. \qquad (B.11) $$

Applying Tadmor's condition using the entropy potential (B.7) and the iden-
tity $[ab] = [a]\bar{b} + \bar{a}[b]$ leads to:

$$[v]_{i+\frac{1}{2}} = g[h]_{i+\frac{1}{2}} \bar{A}_{i+\frac{1}{2}} \bar{u}_{i+\frac{1}{2}} + [\kappa]_{i+\frac{1}{2}} \bar{A}_{i+\frac{1}{2}} \bar{u}_{i+\frac{1}{2}} + [u]_{i+\frac{1}{2}} \bar{A}_{i+\frac{1}{2}} (\bar{u})^2_{i+\frac{1}{2}} + g[u]_{i+\frac{1}{2}} \bar{K}_{i+\frac{1}{2}}, \qquad (B.12)$$

$$= g[h]_{i+\frac{1}{2}} \bar{A}_{i+\frac{1}{2}} \bar{u}_{i+\frac{1}{2}} + g[u]_{i+\frac{1}{2}} \bar{K}_{i+\frac{1}{2}}, \qquad (B.13)$$

$$= g[h]_{i+\frac{1}{2}} \bar{A}_{i+\frac{1}{2}} \bar{u}_{i+\frac{1}{2}} + [v]_{i+\frac{1}{2}} - g\bar{u}_{i+\frac{1}{2}}[K]_{i+\frac{1}{2}}. \qquad (B.14)$$

We assume that $\bar{u}$ is evaluated using the arithmetic average. Thus, the remaining condition to be satisfied is

$$[h]_{i+\frac{1}{2}} \bar{A}_{i+\frac{1}{2}} = [K]_{i+\frac{1}{2}}, \qquad (B.15)$$

which means that $\bar{A}_{i+\frac{1}{2}}$ needs to be an exact mean value linearization of the relation (2.28). A simple definition to obtain an exactly EC numerical flux is to compute $\bar{A}_{i+\frac{1}{2}}$ evaluating exactly the integral

$$\bar{A}_{i+\frac{1}{2}} \doteq \frac{1}{h_{i+\frac{1}{2}}^R - h_{i+\frac{1}{2}}^L} \int_{h_{i+\frac{1}{2}}^L}^{h_{i+\frac{1}{2}}^R} A(h) \, dh. \qquad (B.16)$$

For a symmetric trapezoidal channel for example

$$A(h) = bh + \cot \theta h^2, \qquad (B.17)$$

we then have after few manipulations:

$$\bar{A}_{i+\frac{1}{2}} = \frac{1}{2} (A_{i+\frac{1}{2}}^L + A_{i+\frac{1}{2}}^R) - \frac{1}{6} \cot \theta (h_{i+\frac{1}{2}}^R - h_{i+\frac{1}{2}}^L)^2. \qquad (B.18)$$
## Appendix C. Study of the singularity of the solitary wave solution

The integral (4.13) has a singularity in $h/h_0 = 1+\epsilon$. We studied the function to be integrated near the peak, defining firstly the following quantities:

$$ \zeta \doteq 1 + \epsilon - H, \qquad (C.1) $$

$$ \varrho \doteq \frac{\varphi \psi'}{\psi} (1 + \epsilon) - \varphi'(1 + \epsilon). \qquad (C.2) $$

By expanding to order 1 the function under the square root of (4.13), we obtain:

$$ \varphi(H) - \frac{\varphi(1+\epsilon)}{\psi(1+\epsilon)} \psi(H) \approx \varphi(1+\epsilon) + \varphi'(1+\epsilon) (H - (1+\epsilon)) - \frac{\varphi(1+\epsilon)}{\psi(1+\epsilon)} [\psi(1+\epsilon) + \psi'(1+\epsilon) (H - (1+\epsilon))] \\ = \varrho \zeta. $$

We therefore remove this singularity, and change the variable of integration, so that the integral (4.13) becomes, in the neighbourhood of the peak:

$$ \left( \int_{0}^{1+\epsilon-h/h_0} \left( \frac{1+2\tau_0(1+\epsilon-\zeta)}{\sqrt{\varphi(1+\epsilon-\zeta) - \frac{\varphi(1+\epsilon)}{\psi(1+\epsilon)}\psi(1+\epsilon-\zeta)}} - \frac{1}{\sqrt{\varrho\zeta}} \right) d\zeta \right)_{\zeta=0} \approx \left( \int_{0}^{1+\epsilon-h/h_0} \frac{2\tau_0(1+\epsilon)}{\sqrt{\varrho\zeta}} d\zeta \right)_{\zeta=0} \\ = \left( \frac{4\tau_0(1+\epsilon)}{\sqrt{\varrho}} \left[ \sqrt{\zeta} \right]_{0}^{1+\epsilon-h/h_0} \right)_{\zeta=0} \\ = 0. $$

The function to be integrated diverges in $1+\epsilon$ while we have just shown that its integral converges. The solution that deals with this singularity is therefore:

$$ \frac{x - ct}{h_0} = \mu \sqrt{3(1 + \tau_0) \frac{\gamma^*}{h_0^2}} \left[ \int_{h/h_0}^{1+\epsilon} \left( \frac{1 + 2\tau_0 H}{\sqrt{\varphi(H) - \frac{\varphi(1+\epsilon)}{\psi(1+\epsilon)}\psi(H)}} - \frac{1}{\sqrt{\varrho\zeta}} \right) dH - 2\sqrt{\frac{1 + \epsilon - \frac{h}{h_0}}{\varrho}} \right]. \qquad (C.3) $$
## References

[1] G. Wei, J. T. Kirby, Time-Dependent Numerical Code for Extended Boussinesq Equations, Journal of Waterway, Port, Coastal, and Ocean Engineering 121 (5) (1995) 251-261.
doi:10.1061/(ASCE)0733-950X(1995)121:5(251).
URL http://ascelibrary.org/doi/10.1061/%28ASCE%290733-950X%281995%29121%3A5%28251%29

[2] M. Walkley, M. Berzins, A finite element method for the two-dimensional extended Boussinesq equations, International Journal for Numerical Methods in Fluids 39 (10) (2002) 865-885.
eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/fld.349.
doi:10.1002/fld.349.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.349

[3] C. Eskilsson, S. J. Sherwin, AN hp/SPECTRAL ELEMENT MODEL FOR EFFICIENT LONG-TIME INTEGRATION OF BOUSSINESQ-TYPE EQUATIONS, Coastal Engineering Journal 45 (02) (2003) 295-320, publisher: World Scientific Publishing Co.
doi:10.1142/S0578563403000762.
URL https://www.worldscientific.com/doi/abs/10.1142/S0578563403000762

[4] O. R. Sørensen, H. A. Schäffer, L. S. Sørensen, Boussinesq-type modelling using an unstructured finite element technique, Coastal Engineering 50 (4) (2004) 181-198. doi:10.1016/j.coastaleng.2003.10.005.
URL https://www.sciencedirect.com/science/article/pii/
S0378383903001121

[5] A. P. Engsig-Karup, J. S. Hesthaven, H. B. Bingham, P. A. Mad-
sen, Nodal DG-FEM solution of high-order Boussinesq-type equations,
Journal of Engineering Mathematics 56 (3) (2006) 351–370. doi:
10.1007/s10665-006-9064-z.
URL https://doi.org/10.1007/s10665-006-9064-z

[6] M. Ricchiuto, A. G. Filippini, Upwind residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries, Journal of Computational Physics 271 (2014) 306–341. doi:10.1016/j.jcp.2013.12.048.
URL https://www.sciencedirect.com/science/article/pii/
S0021999113008565

[7] P. Bacigaluppi, M. Ricchiuto, P. Bonneton, Implementation and Evaluation of Breaking Detection Criteria for a Hybrid Boussinesq Model, Water Waves (2020) 207–241 Publisher: Springer. doi:10.1007/s42286-019-00023-8.
URL https://hal.science/hal-02383888

[8] R. Cienfuegos, E. Barthélemy, P. Bonneton, A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: boundary conditions and validation, International Journal for Numerical Methods in Fluids 53 (9) (2007) 1423–1455, eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/fld.1359. doi:10.1002/fld.1359.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.
1359

[9] D. Dutykh, T. Katsaounis, D. Mitsotakis, Finite volume schemes for dis-
persive wave propagation and runup, Journal of Computational Physics
230 (8) (2011) 3035–3061, publisher: Elsevier. doi:10.1016/j.jcp.
2011.01.003.

URL https://hal.science/hal-00472431

[10] H. Ranocha, D. Mitsotakis, D. I. Ketcheson, A Broad Class of
Conservative Numerical Methods for Dispersive Wave Equations,
Communications in Computational Physics 29 (4) (2021) 979–1029,
arXiv:2006.14802 [physics]. doi:10.4208/cicp.OA-2020-0119.

URL http://arxiv.org/abs/2006.14802

[11] P. Bonneton, F. Chazel, D. Lannes, F. Marche, M. Tissier, A splitting approach for the fully nonlinear and weakly dispersive Green–Naghdi model, Journal of Computational Physics 230 (4) (2011) 1479–1498.
doi:10.1016/j.jcp.2010.11.015.

URL https://linkinghub.elsevier.com/retrieve/pii/
S0021999110006236

[12] A. Filippini, M. Kazolea, M. Ricchiuto, A flexible genuinely non-linear approach for nonlinear wave propagation, breaking and run-up, Journal of Computational Physics 310 (2016) 381–417.
doi:10.1016/j.jcp.2016.01.027.

URL https://linkinghub.elsevier.com/retrieve/pii/
S0021999116000358
[13] 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) (2019) 509-531. doi:10.1002/fld.4766.
URL https://onlinelibrary.wiley.com/doi/10.1002/fld.4766

[14] D. Lannes, F. Marche, A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations, Journal of Computational Physics 282 (2015) 238-268, publisher: Elsevier.
doi:10.1016/j.jcp.2014.11.016.
URL https://hal.science/hal-00932858

[15] A. Duran, F. Marche, A discontinuous Galerkin method for a new class of Green-Naghdi equations on simplicial unstructured meshes, Applied Mathematical Modelling 45 (2017) 840-864.
doi:10.1016/j.apm.2017.01.030.
URL https://www.sciencedirect.com/science/article/pii/S0307904X17300367

[16] A. Cauquis, M. Ricchiuto, P. Heinrich, Lax-Wendroff Schemes with Polynomial Extrapolation and Simplified Lax-Wendroff Schemes for Dispersive Waves: A Comparative Study, Water Waves 4 (3) (2022) 345-377. doi:10.1007/s42286-022-00060-w.
URL https://doi.org/10.1007/s42286-022-00060-w

[17] D. Torlo, M. Ricchiuto, Model order reduction strategies for weakly dispersive waves, Mathematics and Computers in Simulation 205 (2023) 997-1028. doi:10.1016/j.matcom.2022.10.034.
URL https://linkinghub.elsevier.com/retrieve/pii/
S0378475422004487

[18] M. Kazolea, A. Filippini, 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 (2023) 102157. doi:10.1016/j.ocemod.2022.102157.
URL https://linkinghub.elsevier.com/retrieve/pii/
S1463500322001718

[19] A. Mazaheri, M. Ricchiuto, H. Nishikawa, A first-order hyperbolic system approach for dispersion, Journal of Computational Physics 321 (2016) 593–605. doi:10.1016/j.jcp.2016.06.001.
URL https://www.sciencedirect.com/science/article/pii/
S0021999116302261

[20] N. Favrie, S. Gavrilyuk, A rapid numerical method for solving Serre-Green-Naghdiequations describing long free surface gravity waves, Nonlinearity Publisher: IOP Publishing (2017).
URL https://hal.science/hal-01358630

[21] S. Tkachenko, S. Gavrilyuk, J. Massoni, Extended Lagrangian approach for the numerical study of multidimensional dispersive waves: applications to the Serre-Green-Naghdi equations, Journal of Computational Physics 477 (2023) 111901, arXiv:2203.10901 [cs, math]. doi:10.1016/j.jcp.2022.111901.
URL http://arxiv.org/abs/2203.10901
[22] S. Gavrilyuk, K.-M. Shyue, Hyperbolic approximation of the BBM equation, Nonlinearity 35 (3) (2022) 1447, publisher: IOP Publishing.
doi:10.1088/1361-6544/ac4c49.
URL https://dx.doi.org/10.1088/1361-6544/ac4c49

[23] S. Busto, M. Dumbser, C. Escalante, N. Favrie, S. Gavrilyuk, On High Order ADER Discontinuous Galerkin Schemes for First Order Hyperbolic Reformulations of Nonlinear Dispersive Systems, Journal of Scientific Computing 87 (2) (2021) 48. doi:10.1007/s10915-021-01429-8.
URL https://doi.org/10.1007/s10915-021-01429-8

[24] C. Escalante, M. Dumbser, M. J. Castro, An efficient hyperbolic relaxation system for dispersive non-hydrostatic water waves and its solution with high order discontinuous Galerkin schemes, Journal of Computational Physics 394 (2019) 385-416. doi:10.1016/j.jcp.2019.05.035.
URL https://www.sciencedirect.com/science/article/pii/S0021999119303730

[25] N. Bonneton, P. Bonneton, J.-P. Parisot, A. Sottolichio, G. Detandt, Ressaut de marée et Mascaret – exemples de la Garonne et de la Seine, Comptes Rendus Geoscience 344 (10) (2012) 508-515.
doi:10.1016/j.crte.2012.09.003.
URL https://linkinghub.elsevier.com/retrieve/pii/S1631071312001551

[26] F. Benet, J. A. Cunge, Analysis Of Experiments On Secondary Undulations Caused By Surge Waves In Trapezoidal Channels, Journal of Hydraulic Research 9 (1) (1971) 11-33. doi:10.1080/00221687109500335.
URL https://www.tandfonline.com/doi/full/10.1080/
00221687109500335

[27] D. Violeau, Contribution to the Theory of Undular Bores. A journey around the Korteweg-de Vries equation, in: IAHR Water Monograph Series, Vol. 1, The International Association for Hydro-Environment Engineering and Research (IAHR), Spain, 2022, pp. 1-96. doi:10.3850/IAHRWaterMonograph001.
URL https://www.iahr.org/library/infor?pid=20693

[28] H. Favre, Étude théorique et expérimentale des ondes de translation dans les canaux découverts, Publications du Laboratoire de recherches hydrauliques annexé à l'École polytechnique fédérale de Zurich, Dunod, 1935.

[29] J. A. Sandover, C. Taylor, Cnoidal waves and bores, La Houille Blanche 48 (3) (1962) 443-465. doi:10.1051/lhb/1962045.
URL https://www.tandfonline.com/doi/full/10.1051/lhb/1962045

[30] A. Treske, Undular bores (favre-waves) in open channels - Experimental studies, Journal of Hydraulic Research 32 (3) (1994) 355-370. doi:10.1080/00221689409498738.
URL https://www.tandfonline.com/doi/full/10.1080/00221689409498738

[31] D. H. Peregrine, Long waves in a uniform channel of arbitrary cross-section, Journal of Fluid Mechanics 32 (2) (1968) 353-365.
doi:10.1017/S0022112068000777.

URL https://www.cambridge.org/core/product/identifier/
S0022112068000777/type/journal_article

[32] J. D. Fenton, Cnoidal waves and bores in uniform channels of arbitrary cross-section, Journal of Fluid Mechanics 58 (3) (1973) 417-434.
doi:10.1017/S0022112073002259.

URL https://www.cambridge.org/core/product/identifier/
S0022112073002259/type/journal_article

[33] N. Shuto, Nonlinear Long Waves in a Channel of Variable Section, Coastal Engineering in Japan 17 (1) (1974) 1-12. doi:10.1080/05785634.1974.11924178.

URL https://www.tandfonline.com/doi/full/10.1080/05785634.1974.11924178

[34] J. W. Miles, On the Korteweg-de Vries equation for a gradually varying channel, Journal of Fluid Mechanics 91 (1) (1979) 181-190, publisher: Cambridge University Press. doi:10.1017/S0022112079000100.

URL https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/on-the-kortewegde-vries-equation-for-a-gradually-varying-channel/C591BDCB524A1A33F597126C0B506436

[35] M. H. Teng, T. Y. Wu, Nonlinear water waves in channels of arbitrary shape, Journal of Fluid Mechanics 242 (1992) 211-233.
doi:10.1017/S0022112092002349.
URL https://www.cambridge.org/core/product/identifier/
S0022112092002349/type/journal_article

[36] M. H. Teng, T. Y. Wu, Evolution of long water waves in vari-
able channels, Journal of Fluid Mechanics 266 (1994) 303–317.
doi:10.1017/S0022112094001011.

URL https://www.cambridge.org/core/product/identifier/
S0022112094001011/type/journal_article

[37] M. H. Teng, T. Y. Wu, Effects of channel cross-sectional geometry on long wave generation and propagation, Physics of Fluids 9 (11) (1997) 3368–3377. doi:10.1063/1.869449.

URL http://aip.scitation.org/doi/10.1063/1.869449

[38] M. H. Teng, Boussinesq solution for solitary waves in uniform channels with sloping side walls, Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science 214 (6) (2000) 781–787. doi:10.1243/0954406001523777.

URL http://journals.sagepub.com/doi/10.1243/
0954406001523777

[39] P. Winckler, P. L.-F. Liu, Long waves in a straight channel with non-uniform cross-section, Journal of Fluid Mechanics 770 (2015) 156–188. doi:10.1017/jfm.2015.147.

URL https://www.cambridge.org/core/product/identifier/
S0022112015001470/type/journal_article

[40] M. A. Debyaoui, M. Ersoy, Generalised Serre–Green–Naghdi equations
for open channel and for natural river hydraulics, Asymptotic Analysis
124 (3-4) (2021) 343-369. doi:10.3233/ASY-201647.

URL https://www.medra.org/servlet/aliasResolver?alias=iospress&doi=10.3233/ASY-201647

[41] D. Violeau, Serre equations in channels and rivers of arbitrary cross section, Journal of Hydraulic EngineeringPublisher: American Society of Civil Engineers (2023).

URL https://hal.science/hal-04050777

[42] U. S. Fjordholm, S. Mishra, E. Tadmor, Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography, Journal of Computational Physics 230 (14) (2011) 5587–5609. doi:10.1016/j.jcp.2011.03.042.

URL https://linkinghub.elsevier.com/retrieve/pii/S0021999111002038

[43] R. Chassagne, A. G. Filippini, M. Ricchiuto, P. Bonneton, Dispersive and dispersive-like bores in channels with sloping banks, Journal of Fluid Mechanics 870 (2019) 595–616. doi:10.1017/jfm.2019.287.

URL https://www.cambridge.org/core/product/identifier/S0022112019002878/type/journal_article

[44] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics: a practical introduction, 3rd Edition, Springer, Dordrecht ; New York, 2009, oCLC: ocn401321914.

[45] L. Brudvik-Lindner, D. Mitsotakis, A. E. Tzavaras, Oscillatory and
regularized shock waves for a dissipative Peregrine-Boussinesq system,
arXiv:2209.10129 [physics] (Mar. 2023). doi:10.48550/arXiv.2209.10129.

URL http://arxiv.org/abs/2209.10129

[46] G. A. El, M. A. Hoefer, Dispersive shock waves and modulation theory, Physica D: Nonlinear Phenomena 333 (2016) 11–65.
doi:10.1016/j.physd.2016.04.006.

URL https://www.sciencedirect.com/science/article/pii/S0167278916301580

[47] P. L. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes, Journal of Computational Physics 43 (2) (1981) 357–372. doi:10.1016/0021-9991(81)90128-5.

URL https://www.sciencedirect.com/science/article/pii/0021999181901285

[48] B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method, Journal of Computational Physics 32 (1) (1979) 101–136. doi:10.1016/0021-9991(79)90145-1.

URL https://www.sciencedirect.com/science/article/pii/0021999179901451

[49] M. E. Vázquez-Cendón, Improved Treatment of Source Terms in Upwind Schemes for the Shallow Water Equations in Channels with Irregular Geometry, Journal of Computational Physics 148 (2) (1999) 497–526. doi:10.1006/jcph.1998.6127.
URL https://linkinghub.elsevier.com/retrieve/pii/
S0021999198961277

[50] E. Tadmor, The Numerical Viscosity of Entropy Stable Schemes for Systems of Conservation Laws. I, Mathematics of Computation 49 (179) (1987) 91–103, publisher: American Mathematical Society. doi:10.2307/2008251.

URL https://www.jstor.org/stable/2008251

[51] D. Mitsotakis, H. Ranocha, D. I. Ketcheson, E. Süli, A conservative fully-discrete numerical method for the regularised shallow water wave equations, SIAM Journal on Scientific Computing 43 (2) (2021) B508–B537, arXiv:2009.09641 [cs, math]. doi:10.1137/20M1364606.

URL http://arxiv.org/abs/2009.09641

[52] L. Thomas, Elliptic Problems in Linear Difference Equations over a Network (1949).

[53] F. Bouchut, An introduction to finite volume methods for hyperbolic conservation laws, ESAIM: Proceedings 15 (2005) 1–17, publisher: EDP Sciences. doi:10.1051/proc:2005020.

URL https://www.esaim.proc.org/articles/proc/abs/2005/02/bouchut/bouchut.html

[54] S. Gottlieb, C.-W. Shu, E. Tadmor, Strong Stability-Preserving High-Order Time Discretization Methods, SIAM Review 43 (1) (2001) 89–112, publisher: Society for Industrial and Applied Mathematics.

URL https://www.jstor.org/stable/3649684
[55] M.-O. Bristeau, B. Coussin, Boundary Conditions for the Shallow Water Equations solved by Kinetic Schemes, report, INRIA (2001).
URL https://inria.hal.science/inria-00072305

[56] K. Salari, P. Knupp, Code Verification by the Method of Manufactured Solutions, Tech. Rep. SAND2000-1444, Sandia National Lab. (SNL-NM), Albuquerque, NM (United States); Sandia National Lab. (SNL-CA), Livermore, CA (United States) (Jun. 2000). doi:10.2172/759450.
URL https://www.osti.gov/biblio/759450

[57] P. J. Roache, Code Verification by the Method of Manufactured Solutions, Journal of Fluids Engineering 124 (1) (2001) 4-10. doi: 10.1115/1.1436090.
URL https://doi.org/10.1115/1.1436090

[58] F. Serre, Contribution à l'étude des écoulements permanents et variables dans les canaux, La Houille Blanche 39 (6) (1953) 830-872.
doi:10.1051/lhb/1953058.
URL https://www.tandfonline.com/doi/full/10.1051/lhb/1953058

[59] C. C. Mei, The applied dynamics of ocean surface waves, Ocean Engineering 11 (3) (1984) 321. doi:10.1016/0029-8018(84)90033-7.
URL https://www.sciencedirect.com/science/article/pii/0029801884900337