# Exploring different possibilities for second-order well-balanced Lagrange-projection numerical schemes applied to shallow water Exner equations

Christophe Chalons, Alessia Del Grosso

## ► To cite this version:

Christophe Chalons, Alessia Del Grosso. Exploring different possibilities for second-order well-balanced Lagrange-projection numerical schemes applied to shallow water Exner equations. International Journal for Numerical Methods in Fluids, 2022, 94 (6), pp.505-535. 10.1002/fld.5064. hal-03993015

HAL Id: hal-03993015

https://hal.science/hal-03993015

Submitted on 28 Feb 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# Exploring different possibilities for second-order well-balanced Lagrange-projection numerical schemes applied to shallow water Exner equations

C. Chalons*and A. Del Grosso†

**Abstract.** This work is devoted to the numerical approximation of the shallow water Exner system. We investigate three different numerical strategies to discretize the Exner equation which expresses the evolution in time of the bed sediment. The numerical schemes are all based on the Lagrange-projection formalism which consists in splitting the mathematical model into the acoustic and transport system. In particular, the Exner equation is taken into account either in the acoustic and transport steps, or only at the acoustic or transport level. The methods and their second-order extensions are designed in such a way to satisfy the well-balanced property, namely the "lake at rest" and the "constant bed slope" steady states. Numerical evidences are given to validate the numerical schemes.

## 1 Introduction and governing equations

This work considers several second-order and well-balanced Lagrange-projection schemes applied to the shallow water system with moving topography. As it will be seen, Lagrange-projection schemes consist in splitting the acoustic and transport waves of the model in two different systems (and steps), expedient that can be very useful in practice, for instance in subsonic regimes. This kind of decomposition can also be interpreted as a Lagrange-projection one, in the sense that the considered system is first written in Lagrangian coordinates and solved. Then, the Lagrangian solution is projected again into Eulerian coordinates. Here we do not give further details and refer the reader for instance to [15, 14, 5, 11, 16] for first-order Lagrange-projection scheme and to [37, 21, 9] for methods of higher order of accuracy. Our main objective here is to focus on the topography discretization and to compare several natural approaches in which the bed level is taken into account either partially in both steps, or entirely in one of the two steps. Notice that, in the previous work [13], it has been described a 1D-2D second-order well-balanced Lagrange-projection scheme for the shallow water Exner system with bed level discretization completely taken into account in the transport step, resembling a splitting method. As such, the present work is considered as its natural extension.

<footnote>*Université Paris-Saclay, UVSQ, CNRS, Laboratoire de Mathématiques de Versailles, 78000, Versailles, France. E-mail: christophe.chalons@uvsq.fr</footnote>

<footnote>†Université Paris-Saclay, UVSQ, CNRS, Laboratoire de Mathématiques de Versailles, 78000, Versailles, France. E-mail: alessia.del-grosso2@uvsq.fr</footnote>
We underline that often the Saint-Venant equations with moving topography have been numerically solved by means of splitting methods meaning that the hydrodynamic and morphologic components are separately solved. Indeed, such methods are easier to implement even if they are known to produce spurious oscillations in the numerical results. The latter are mainly due to reasons related to the hyperbolic properties of the coupled systems, for which we directly refer to [20] for further details. On the other hand, also coupled Riemann solvers have been proposed for the 1D-2D shallow water Exner system able to remove (at least in part) the above mentioned oscillations, see for instance [39] where a Roe-type first-order scheme was used or [7] in which a path-conservative Roe method and its high-order extension together with flux limiters has been presented. Another interesting reference is [33], where two schemes based on the Roe approach have been described. In particular, the first one is based on a fully coupled approach while the second one consists in a decoupled strategy where the oscillations are in part stabilized by controlling the stability region. Without being exhaustive, see for instance [22, 32, 40] for further coupled schemes applied to the shallow water Exner system.

Let us now present the model we are interested in, namely the shallow water Exner system. The first two equations have been extensively used to describe the evolution in time of fluid flows for instance in rivers or coastal areas. Whereas, the third equation simulates the bedload sediment transport due to the mechanical action of the fluid. As such, the system reads

$$ \left\{ \begin{array}{l} \partial_t h + \partial_x (hu) = 0 \\ \partial_t (hu) + \partial_x (hu^2 + p) + gh\partial_x z = -ghS_f \\ \partial_t z + \zeta\partial_x q_b = 0, \end{array} \right. \qquad (1.1) $$

where $h(x, t) > 0$ is the water depth, $u(x, t)$ the averaged velocity, $z(x, t)$ the bed level and, in particular, $H = h + z$ is the free surface elevation. Then, $p = \frac{gh^2}{2}$ is the pressure term with $g$ the gravitational acceleration, $q_b = q_b(h, u)$ is the solid transport discharge and $\zeta$ a parameter such that $\zeta = \frac{1}{1-\rho_0}$ with $\rho_0$ the porosity of the sediment layer. Moreover, $-ghS_f$ represents the Manning friction term where $S_f = \frac{\mu_f^2|u|u}{R_h^{4/3}}$ with $\mu_f$ Manning roughness coefficient and $R_h = \frac{Lh}{L+2h}$ hydraulic radius, where $L$ is the length of the channel. Observe that this definition of $R_h$ is valid for rectangular channel. Finally, $t > 0$ represents the time and $x$ the axial coordinate. In compact form, we have

$$ \partial_t \mathbf{Q} + \partial_x \mathbf{F}(\mathbf{Q}) + \mathbf{A}(\mathbf{Q})\partial_x \mathbf{Q} = -ghS_f \mathbf{E}_2 $$

where $\mathbf{E}_2 = (0, 1, 0)^t$,

$$ \mathbf{Q} = \begin{pmatrix} h \\ hu \\ z \end{pmatrix}, \quad \mathbf{F}(\mathbf{Q}) = \begin{pmatrix} hu \\ hu^2 + p \\ \zeta q_b \end{pmatrix}, \quad \mathbf{A}(\mathbf{Q}) = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & gh \\ 0 & 0 & 0 \end{pmatrix}. $$

For more details about shallow-water equations with and without moving topography, we refer for instance to [1, 6, 7, 4] and [42, 2, 8, 11, 38] [34, 32]. Let us know focus on the solid transport discharge $q_b$, which can be formulated in different ways depending on the characteristics of the sediment and the flow, see for instance [27]. One frequently used formulation is the well-known Grass model, which expresses the instantaneous sediment transport as a power law of the averaged velocity $u$, namely

$$ q_b = A_g u |u|^{m_g-1}, \quad 1 \le m_g \le 4, \qquad (1.2) $$
where the parameter $A_g \in [0, 1]$ is computed using empirical relationships based on the local properties (grain size, cinematic viscosity..). Moreover, it expresses how strong is the interaction between the flow and the sediment. For instance, the interaction is considered weak if $A_g$ is of order $10^{-3}$ or smaller, while for values of order $10^{-1}$ the flow is said to be highly erosive, see for instance [7, 33, 39, 31, 3] for further details. In practice, we will set $m_g = 3$ and we briefly recall that the Grass closure relation leads to a strictly hyperbolic system with all real eigenvalues, see again [7]. It is important to underline that, when using the Grass formulation, we are implicitly supposing that the bed sediments start moving as soon as the velocity of the water is different from zero. However, in other formulations it is usually assumed that a critical value has to be overcome. Let us see the details considering for instance the Meyer-Peter&Müller (MPM) formulation. The latter is given by

$$ q_b = 8Qsgn(u)(\theta^* - \theta_c^*)_+^{\frac{3}{2}} \quad \text{with} \quad \theta^* = \frac{u_*^2}{sgd} \quad \text{and} \quad u_*^2 = \frac{g\mu_f^2 u^2}{h^{\frac{1}{3}}} \quad (1.3) $$

with $Q = d\sqrt{gsd}$ the characteristic discharge where $s$ is the relative density and $d$ the sediment diameter. Then, $\theta^*$ represents the non-dimensional shear stress and determines the movement of the sediments. Indeed, only if it is bigger of the critical stress value $\theta_c^*$, $q_b$ is different from zero. Here we do not present further formulations but many have been proposed, see again the previous references. We also highlight that each formulation has usually its own range of applications which depends on the characteristics of the flow and sediments. For instance the MPM formula is only used for weakly erosive flow, see directly [7] for details about the range of values for the parameters. Finally, let us specify that, depending on the particular form of $q_b$, the convective part of system (1.1) could be strictly hyperbolic or not.

Last but not least, we are interested in numerical schemes able to preserve the stationary solutions of the system. This property is in general not trivial to satisfy but, at the same time, critical if we want to obtain accurate numerical methods which do not produce spurious oscillations in their results when near to stationary solutions. If such a property is met, the numerical scheme is said to be well-balanced, which, in our specific case, means the preservations of the following steady states,

$$ q = hu = \text{constant} = q_0, \quad \partial_x \left( \frac{q_0^2}{2h^2} + g(h+z) \right) + gS_f = 0 \quad \text{and} \quad q_b = \text{constant}. $$

Here, we are only interested in preserving two particular steady states, namely the so-called "lake at rest" equilibrium with zero-velocity

$$ u = 0 \quad \text{and} \quad h + z = \text{constant}, \quad (1.4) $$

and the "constant bed slope" equilibrium

$$ \partial_x h = \partial_x u = \partial_{xx} z = 0, \quad \text{and} \quad \partial_x z + S_f = 0, \quad (1.5) $$

Let us observe that when using the Grass formulation, these are the only two possible steady states. Being the well-balancedness of the scheme a crucial property, many studies have been done in this sense, here we refer for instance to [2, 3, 8, 30] and the references herein. See also [11, 16, 37] for well-balanced methods in the Lagrange-projection formalism.
*Outline of the paper.* We now give the structure of the paper. In the next section we briefly present the Lagrange-projection splitting strategy in both Eulerian and Lagrangian variables for the usual shallow water system, thus for the evolution equations of $h$ and $hu$. Then, in section 3 we explain three different strategies to treat the topography equation. In particular, details for the approximate Riemann solver for the acoustic systems are given. Subsequently, we present the numerical schemes both at first and second-order of accuracy in sections 4 and 5 respectively. Section 6 is devoted to the description of the well-balanced property for each scheme. Finally, we exhibit several numerical evidences to validate our numerical schemes in section 7, while concluding remarks are given in section 8.

## 2 Operator splitting for the shallow water system

This section focuses only on the first two equations of system (1.1), namely the updating equations for the water height $h$ and discharge $hu$, whereas details for the bed level evolution will be given later. Thus, here we explain the decomposition which entails the splitting of the Saint-Venant system into two different ones, the so-called acoustic and transport systems. The former accounts for the acoustic phenomena and topography variations, while the latter focuses on the transport effects. Note that here we neglect the friction term; its contribution will be included directly at the end of the numerical methods, see section 4.

Then, we observe that the first two equations of the model can be reformulated as

$$ egin{cases} \partial_t h + h \partial_x u + u \partial_x h = 0 \\ \partial_t (hu) + h u \partial_x u + u \partial_x (hu) + \partial_x \left(\frac{gh^2}{2}\right) = -gh \partial_x z, \end{cases} $$

where we used the chain rule for space derivatives. Therefore, the so-called acoustic and transport systems are respectively given by

$$ egin{cases} \partial_t h + h \partial_x u = 0 \\ \partial_t (hu) + h u \partial_x u + \partial_x \left(\frac{gh^2}{2}\right) = -gh \partial_x z, \end{cases} \quad (2.1) $$

and

$$ egin{cases} \partial_t h + u \partial_x h = 0 \\ \partial_t (hu) + u \partial_x (hu) = 0, \end{cases} $$

where the latter can be reinterpreted as

$$ \partial_t X + u \partial_x X = 0 $$

with $X$ either $X = h$ or $X = hu$. We also observe that system (2.1) can be expressed as

$$ egin{cases} \partial_t \tau - \partial_m u = 0 \\ \partial_t u + \partial_m p = - \frac{g}{\tau} \partial_m z \end{cases} $$

where $\tau = \frac{1}{h}$ and the mass variable $m$ is such that $\frac{1}{h} \partial_x = \partial_m$. See for instance [11, 16] for additional details about this decomposition applied to the shallow water system. It is then clear that the numerical
method would sum up to first solve the acoustic system and, after that, the transport equations with the acoustic solution as initial data. However, while it is well known how to decompose the shallow water system, this is not true when considering the Exner equation. Indeed, one could easily imagine at least three possibilities for numerical treatment. The first one would split it inside both steps as we just did for h and hu, the second one would account for it at the acoustic level, and the last one directly inside the transport step. The aim of the present contribution is indeed to compare these three approaches, both at first and second order accuracy.

That being said and before going into further details, it is convenient to first introduce the Lagrangian coordinates. We first define the fluid particle $\xi$ and the characteristic curves

$$ \begin{cases} \frac{\partial x}{\partial t}(\xi, t) = u(x(\xi, t), t) \\ x(\xi, 0) = \xi \end{cases} $$

which define the trajectory: $t \to x(\xi, t)$, of $\xi$ as the time goes on. Therefore, any function: $(x, t) \to \varphi(x, t)$ in Eulerian coordinates can be written in Lagrangian coordinates,

$$ \bar{\varphi}(\xi, t) = \varphi(x(\xi, t), t). $$

Introducing now the volume ratio

$$ L(\xi, t) = \frac{\partial x}{\partial \xi}(\xi, t) $$

such that

$$ \begin{cases} \frac{\partial L}{\partial t}(\xi, t) = \partial_{\xi} u(x(\xi, t), t) \\ L(\xi, 0) = 1, \end{cases} \qquad (2.2) $$

it clearly follows

$$ \partial_t L(\xi, t) = \partial_{\xi} u(x(\xi, t), t) = \partial_{\xi} \bar{u}(\xi, t), $$

and thus

$$ \partial_{\xi} \bar{\varphi}(\xi, t) = L(\xi, t) \partial_{x} \varphi(x, t) \quad \text{and} \quad \partial_{t} \bar{\varphi}(\xi, t) = \partial_{t} \varphi(x, t) + u(x, t) \partial_{x} \varphi(x, t). $$

Focusing first on the governing equations for h and hu, it is easy to show that their counterpart in Lagrangian coordinates reads

$$ \begin{cases} \partial_t (L \bar{h}) = 0 \\ \partial_t (L \bar{h} u) + \partial_\xi \bar{p} = -g \bar{h} \partial_\xi \bar{z}. \end{cases} \qquad (2.3) $$

In the following sections, we shall sometimes omit the bar over the Lagrangian functions to avoid cumbersome notations. Hence, the two-steps (acoustic and transport) numerical method would now consist of solving the Lagrangian system (2.3) and then projecting the solution into Eulerian coordinates. For further details about this decomposition, we refer the reader to [37] and the references therein. Let us now discuss the three different strategies proposed for the Exner equations.
## 3 Treatment of the Exner equation

The treatment of the Exner equation is an important issue due to the complexity of the fully coupled system. It is known that a fully decoupled scheme, which consists in updating the topography independently from the first two equations of the model generally produces spurious oscillations inside the numerical solutions, see for instance [20]. However and in order to avoid this problem, a weak coupling of the equations can lead to satisfying numerical results, see [3]. In this work, we mainly focus on weakly coupled numerical schemes. In brief, we are going to take into account the flow and sediment interactions in three ways:

1. The usual acoustic-transport splitting is considered for the topography equation and therefore the bed level is taken into account in both steps;

2. The topography is accounted for only at the level of the Lagrangian step;

3. The topography is updated only in the transport step. This approach has already been proposed in [13] and it resembles the usual splitting (and therefore decoupled) method for shallow water Exner system, but as we will see, spurious oscillations are often not observed here thanks to the Lagrange-projection strategy.

### 3.1 Update the bed level in both steps

The first strategy consists in splitting the bed level evolution equation

$$ \partial_t z + \zeta \partial_x q_b = 0 $$

into the following two equations, namely

$$ \partial_t z - u \partial_x z + \zeta \partial_x q_b = 0 $$

and

$$ \partial_t z + u \partial_x z = 0, $$

so that the complete acoustic system is now given by

$$ \left\{ \begin{array}{l} \partial_t h + h \partial_x u = 0 \\ \partial_t (hu) + hu \partial_x u + \partial_x \left( \frac{gh^2}{2} \right) + gh \partial_x z = 0 \\ \partial_t z - u \partial_x z + \zeta \partial_x q_b = 0 \end{array} \right. \qquad (3.1) $$

while the transport system is simply formulated as

$$ \partial_t X + u \partial_x X = 0 \qquad (3.2) $$

where now we do not only have $X = h, hu$ but also $X = z$. It is thus clear that, with this approach, the evolution of the bed level is taken into account in both the acoustic and transport steps. Then, note
that the acoustic system (3.1) can be easily reformulated as in the following

$$ egin{cases} 
\partial_t \tau - \partial_m u = 0 \\
\partial_t u + \partial_m p + \frac{g}{\tau} \partial_m z = 0 \\
\partial_t z - \frac{u}{\tau} \partial_m z + \zeta \frac{1}{\tau} \partial_m q_b = 0
\end{cases} 
\quad (3.3) $$

exploiting once again the notation $\tau = 1/h$ and $\tau \partial_x = \partial_m$. Alternatively, the shallow water Exner system (1.1) in Lagrangian coordinates reads

$$ egin{cases} 
\partial_t(L\bar{h}) = 0 \\
\partial_t(L\bar{h}u) + \partial_\xi \bar{p} + g\bar{h}\partial_\xi \bar{z} = 0 \\
\partial_t(L\bar{z}) - \partial_\xi(\bar{z}u) + \zeta\partial_\xi \bar{q}_b = 0.
\end{cases} 
\quad (3.4) $$

### 3.1.1 Approximate Riemann solver

In order to approximate the solutions of system (3.4) using a Godunov-type method, we define in this section an approximate Riemann solution of equations (3.3) associated with the initial data

$$ (\tau, u, z)^T(x, t=0) = \begin{cases} (\tau_L, u_L, z_L)^T & \text{if } x < 0 \\ (\tau_R, u_R, z_R)^T & \text{if } x > 0. \end{cases} $$

Here, the idea is to base the approximate Riemann solver on a relaxation formulation, see for instance [10, 16] for more details. Thus, as a starting point, we considered the following relaxation system of the whole (not only acoustic) Saint-Venant-Exner system (1.1),

$$ egin{cases} 
\partial_t h + \partial_x(hu) = 0 \\
\partial_t(hu) + \partial_x(hu^2 + \Pi) + gh\partial_x z = 0 \\
\partial_t z + \partial_x \Omega = 0 \\
\partial_t \Pi + u\partial_x \Pi + a^2\partial_x u = 0 \\
\partial_t \Omega + 2u\partial_x \Omega + (b^2 - \frac{h^2}{u^2})\partial_x z = 0,
\end{cases} 
\quad (3.5) $$

which was proposed in [1]. Observe that we introduced two new variables, $\Pi$ and $\Omega$ which can be respectively interpreted as a linearization of $p$ and $\zeta q_b$. Then, $a$ and $b$ are two constant parameters. However, in order to build a relaxation system for the acoustic system and, at the same time, to define a *well-balanced* approximate Riemann solver, it was convenient to modify system (3.5) in the following way

$$ egin{cases} 
\partial_t \tau - \partial_m u = 0 \\
\partial_t u + \partial_m \Pi = -\frac{g}{\tau} \partial_m z \\
\partial_t z - \frac{u}{\tau} \partial_m z + \frac{1}{\tau} \partial_m \Omega = 0 \\
\partial_t \Pi + a^2 \partial_m u = 0 \\
\partial_t \Omega + \frac{u}{\tau} \partial_m \Omega + u^2 \left(b^2 \tau - \frac{1}{\tau}\right) \partial_m z = 0.
\end{cases} 
\quad (3.6) $$
Thus the idea would be to take the exact Riemann solution of the latter relaxation system (3.6) as the approximate solution of either equations (3.3) or system (3.4). We also underline that now the initial data are

$$ (\tau, u, z, \Pi, \Omega)^T(x, t=0) = \begin{cases} (\tau_L, u_L, z_L, \Pi_L, \Omega_L)^T & \text{if } x < 0 \\ (\tau_R, u_R, z_R, \Pi_R, \Omega_R)^T & \text{if } x > 0 \end{cases} $$

and they are taken at equilibrium, that is to say such that

$$ \Pi_{L,R} = \frac{g}{2} h_{L,R}^2 \quad \text{and} \quad \Omega_{L,R} = \zeta(q_b)_{L,R}. $$

Furthermore, the parameters $a$ and $b$ are chosen in such a way to ensure the stability of the relaxation system. In particular, we ask for the validity of the so-called sub-characteristic condition

$$ a \geq h \sqrt{\partial_h p} \quad \text{and} \quad u^2 b^2 \geq (hu)^2 + gh^2 \zeta \partial_u q_b, \quad b > 0, \qquad (3.7) $$

which is clearly different from the one used for the original relaxation system (3.5). Notice that, in the relaxation system (3.6), which can be written in compact form as

$$ \partial_t \mathbf{U} + \mathbf{B}(\mathbf{U}) \partial_m \mathbf{U} = \mathbf{S}(\mathbf{U}) $$

with

$$ \mathbf{U} = \begin{pmatrix} \tau \\ u \\ z \\ \Pi \\ \Omega \end{pmatrix}, \quad \mathbf{B}(\mathbf{U}) = \begin{pmatrix} 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & -\frac{u}{\tau} & 0 & \frac{1}{\tau} \\ 0 & a^2 & 0 & 0 & 0 \\ 0 & 0 & u^2(b^2\tau - \frac{1}{\tau}) & 0 & \frac{u}{\tau} \end{pmatrix}, \quad \text{and} \quad \mathbf{S}(\mathbf{U}) = \begin{pmatrix} 0 \\ \tilde{s} \\ 0 \\ 0 \\ 0 \end{pmatrix}, $$

the topography term

$$ \tilde{s} = -\frac{g}{\tau} \partial_m z $$

is taken into account as a source term and not included in the convective matrix $\mathbf{B}$. This is also a critical point in order to be able to define a *well-balanced* approximate Riemann solver, even if it is also the reason why the resulting numerical method will only be weakly coupled and not fully coupled. With this in mind, it is easy to show that the matrix $\mathbf{B}$ has five real eigenvalues given by $\lambda_0 = 0, \lambda_a^\pm = \pm a, \lambda_b^\pm = \pm |u|b$, and that the associated characteristic fields are all linearly degenerate. This property is well-known to make the resolution of the Riemann problem straightforward using the continuity of the Riemann invariants across each wave, see for instance [1] for more details. Notice, however, that in this case the eigenvalues are not ordered a priori, so that at a continuous level there exist two different cases depending on whether $a < |u|b$ or not (recall that $a$ and $b$ are positive). In practice, we will distinguish between the following two cases $a < |u_L|b, a < |u_R|b$ and its negation.
### 3.1.2 The case $a < |u_L|b$ and $a < |u_R|b$

In this first case, easy calculations show that the solution of the Riemann problem reads

$$ \hat{U}(\frac{m}{t}; \mathbf{U}_L, \mathbf{U}_R) = \begin{cases} \mathbf{U}_L & \text{if } \frac{m}{t} < \lambda_b^- = -|u_L|b \\ \mathbf{U}_{b,L}^* & \text{if } -|u_L|b < \frac{m}{t} < \lambda_a^- = -a \\ \mathbf{U}_{a,L}^* & \text{if } -a < \frac{m}{t} < \lambda_0 = 0 \\ \mathbf{U}_{a,R}^* & \text{if } 0 < \frac{m}{t} < \lambda_a^+ = a \\ \mathbf{U}_{b,R}^* & \text{if } a < \frac{m}{t} < \lambda_b^+ = |u_R|b \\ \mathbf{U}_R & \text{if } \frac{m}{t} > |u_R|b. \end{cases} $$

with

$$ \mathbf{U}_{b,L}^* = \begin{pmatrix} \tau_L \\ u_L \\ z^* \\ \Pi_L \\ \Omega^* \end{pmatrix}, \quad \mathbf{U}_{a,L}^* = \begin{pmatrix} \tau_L^* \\ u^* \\ z^* \\ \Pi_L^* \\ \Omega^* \end{pmatrix}, \quad \mathbf{U}_{a,R}^* = \begin{pmatrix} \tau_R^* \\ u^* \\ z^* \\ \Pi_R^* \\ \Omega^* \end{pmatrix}, \quad \text{and} \quad \mathbf{U}_{b,R}^* = \begin{pmatrix} \tau_R \\ u_R \\ z^* \\ \Pi_R \\ \Omega^* \end{pmatrix}, \quad (3.8) $$

and on the first hand

$$ \left\{ \begin{aligned} \tau_L^* &= \tau_L + \frac{1}{a}(u^* - u_L) \\ \tau_R^* &= \tau_R - \frac{1}{a}(u^* - u_R) \\ u^* &= \frac{1}{2}(u_L + u_R) - \frac{1}{2a}(\Pi_R - \Pi_L) - \frac{M}{2a} \\ \Pi_{*L} &= \frac{1}{2}(\Pi_L + \Pi_R) - \frac{a}{2}(u_R - u_L) + \frac{M}{2} \\ \Pi_{*R} &= \frac{1}{2}(\Pi_L + \Pi_R) - \frac{a}{2}(u_R - u_L) - \frac{M}{2} \end{aligned} \right. \quad (3.9) $$

where the discretization $M$ of the topography source term reads

$$ M = \frac{g}{2} \left( \frac{1}{\tau_L} + \frac{1}{\tau_R} \right) (z_R - z_L), \quad (3.10) $$

and on the other hand

$$ \left\{ \begin{aligned} z^* &= \frac{|u_R|(\operatorname{sign}(u_R) + b\tau_R)z_R - |u_L|(\operatorname{sign}(u_L) - b\tau_L)z_L}{|u_R|(\operatorname{sign}(u_R) + b\tau_R) - |u_L|(\operatorname{sign}(u_L) - b\tau_L)} - \frac{\Omega_R - \Omega_L}{|u_R|(\operatorname{sign}(u_R) + b\tau_R) - |u_L|(\operatorname{sign}(u_L) - b\tau_L)} \\ \Omega^* &= \frac{\Omega_R + \Omega_L}{2} + \frac{|u_R|}{2}(\operatorname{sign}(u_R) + b\tau_R)(z^* - z_R) + \frac{|u_L|}{2}(\operatorname{sign}(u_L) - b\tau_L)(z^* - z_L). \end{aligned} \right. \quad (3.11) $$
### 3.1.3 The case $a \ge |u_L|b$ or $a \ge |u_R|b$

In this case, the structure of the Riemann solution is the same but the waves are expected to be ordered in a different way. More precisely, assuming that $a > |u^*|b$, we have now

$$ \hat{\mathbf{U}}\left(\frac{m}{t}; \mathbf{U}_L, \mathbf{U}_R\right) = \begin{cases} \mathbf{U}_L & \text{if } \frac{m}{t} < -a \\ \mathbf{U}_{a,L}^* & \text{if } -a < \frac{m}{t} < -u^*b \\ \mathbf{U}_{b,L}^* & \text{if } -u^*b < \frac{m}{t} < 0 \\ \mathbf{U}_{b,R}^* & \text{if } 0 < \frac{m}{t} < u^*b \\ \mathbf{U}_{a,R}^* & \text{if } u^*b < \frac{m}{t} < a \\ \mathbf{U}_R & \text{if } \frac{m}{t} > a. \end{cases} $$

where

$$ \mathbf{U}_{a,L}^* = \begin{pmatrix} \tau_L^* \\ u^* \\ z_L \\ \Pi_L^* \\ \Omega_L^* \end{pmatrix}, \quad \mathbf{U}_{b,L}^* = \begin{pmatrix} \tau_L^* \\ u^* \\ z^* \\ \Pi_L^* \\ \Omega^* \end{pmatrix}, \quad \mathbf{U}_{b,R}^* = \begin{pmatrix} \tau_R^* \\ u^* \\ z^* \\ \Pi_R^* \\ \Omega^* \end{pmatrix}, \quad \text{and} \quad \mathbf{U}_{a,R}^* = \begin{pmatrix} \tau_R^* \\ u^* \\ z_R \\ \Pi_R^* \\ \Omega_R^* \end{pmatrix}. \quad (3.12) $$

Then, (3.9) and (3.10) are still valid, but (3.11) is replaced by

$$ \left\{ \begin{aligned} (z|u|)^* &= |u^*| \frac{(\operatorname{sign}(u^*) + b\tau_R^*)z_L - (\operatorname{sign}(u^*) - b\tau_L^*)z_L}{b(\tau_R^* + \tau_L^*)} - \frac{\Omega_R - \Omega_L}{b(\tau_R^* + \tau_L^*)} \\ \Omega^* &= \frac{\Omega_R + \Omega_L}{2} + \frac{1}{2}((\operatorname{sign}(u^*) + b\tau_R^*)(z|u|)^* - |u^*|z_R) + (\operatorname{sign}(u^*) - b\tau_L^*)(z|u|)^* - |u^*|z_L). \end{aligned} \right. \quad (3.13) $$

Let us notice that assuming $a \ge |u_L|b$ or $a \ge |u_R|b$ does not necessarily imply that $a > |u^*|b$. In practice, we proceed as in the following. First, we define $a$ and $b$ by

$$ a = \max(\Delta x, h_L c_L, h_R c_R) \quad \text{and} \quad b = \max\left(\Delta x, \sqrt{h_L^2 + g \frac{h_L^2}{u_L^2} \zeta \partial_u(q_b)_L}, \sqrt{h_R^2 + g \frac{h_R^2}{u_R^2} \zeta \partial_u(q_b)_R}\right), \quad (3.14) $$

$c = \sqrt{\partial_h p}$, as a natural approximation of (3.7). Then, if $a$ and $b$ are such that $a \ge |u_L|b$ or $a \ge |u_R|b$ but $a \le |u^*|b$, we increase the value of $a$ and redefine it as $a = (1 + \epsilon)|u^*|b$ (with typically $\epsilon = 0.01$). We underline that, once we have redefined $a$, we have to recompute the value of $u^*$, and more generally the quantities in (3.8) and (3.9). In practice, this iterative process converges in one or two iterations. An easier (and more diffusive) option could be to define $a$ and $b$ such that $a$ is automatically smaller than $ub$, in this way we only had to use the star values in section 3.1.2 and no further details would be required. For instance, this is what have been done in [1] for numerically solving system (3.5). However, we will see later in section 6.1 that, to distinguish between the two cases, allows us to obtain a well-balanced numerical scheme.

Notice also that unlike (3.11), we define $(z|u|)^*$ instead of $z^*$ in (3.13) in order to avoid any possible ambiguity related to the value of $u^*$ which can be zero. As we will see below, this is sufficient as in the resulting scheme, the update formula of the topography $z^{n+1-}$ only requires $(uz)^*$ and not $z^*$.
At last, observe that for $u = 0$, the relaxation system (3.6) is not strictly hyperbolic, unlike the La-grangian system (3.3). However, the values in the star region result to be well-defined anyway. Indeed, if $u_L = u_R = u^* = 0$, then (3.13) gives $(uz)^* = 0$ and $\Omega^* = 0$ as well, as $u_L = u_R = 0$ implies $\Omega_L = \Omega_R = 0$.

## 3.2 Update the bed level in the acoustic step

The second strategy takes into account the Exner equation only in the acoustic step. As a consequence, we will have a different acoustic system and therefore a different approximate Riemann solver. Let us give more details.

The acoustic system is now given by

$$ egin{cases} \partial_t h + h \partial_x u = 0 \\ \partial_t (hu) + h u \partial_x u + \partial_x \left(\frac{gh^2}{2}\right) + gh \partial_x z = 0 \\ \partial_t z + \zeta \partial_x q_b = 0 \end{cases} \qquad (3.15) $$

while the transport system is simply formulated as $\partial_t X + u \partial_x X = 0$ with $X = h, hu$, together with $\partial_t z = 0$. Exploiting the variables $\tau = 1/h$ and $\tau \partial_x = \partial_m$, the acoustic system (3.15) also reads

$$ egin{cases} \partial_t \tau - \partial_m u = 0 \\ \partial_t u + \partial_m p + \frac{g}{\tau} \partial_m z = 0 \\ \partial_t z + \zeta \frac{1}{\tau} \partial_m q_b = 0. \end{cases} \qquad (3.16) $$

### 3.2.1 Relaxation system and approximate Riemann Solver

We proceed as before to approximate the solutions of (3.16) and to define an approximate Riemann solution based on a relaxation system. The latter is now defined as

$$ egin{cases} \partial_t \tau - \partial_m u = 0 \\ \partial_t u + \partial_m \Pi = -\frac{g}{\tau} \partial_m z \\ \partial_t z + \frac{1}{\tau} \partial_m \Omega = 0 \\ \partial_t \Pi + a^2 \partial_m u = 0 \\ \partial_t \Omega + u^2 b^2 \tau \partial_m z = 0, \end{cases} $$

with once again

$$ (\tau, u, z, \Pi, \Omega)^T(x, t=0) = \begin{cases} (\tau_L, u_L, z_L, \Pi_L, \Omega_L)^T & \text{if } x < 0 \\ (\tau_R, u_R, z_R, \Pi_R, \Omega_R)^T & \text{if } x > 0 \end{cases} $$

the initial data taken at equilibrium, that is to say such that

$$ \Pi_{L,R} = \frac{g}{2} h_{L,R}^2 \quad \text{and} \quad \Omega_{L,R} = \zeta(q_b)_{L,R}. $$
As far as the parameters a and b are concerned, we now ask for

$$a \geq h\sqrt{\partial_{h}p} \quad \text{and} \quad u^{2}b^{2} \geq gh^{2}\zeta\partial_{u}q_{b}, \quad b > 0. \qquad (3.17)$$

Considering again the topography term in the second equation as a source term, it is easy to show that the characteristic velocities of the model are still given by $\lambda_0 = 0, \lambda_a^{\pm} = \pm a, \lambda_b^{\pm} = \pm|u|b$, and that the associated characteristic fields are all linearly degenerate again. Let us now discuss the associated Riemann solution.

### 3.2.2 Modification of the well-balanced approximate Riemann solver

The procedure is exactly the same as the one presented in section 3.1.1, so that we do not give other details than the differences. In particular, we have to distinguish between the case $a < |u_L|b$ and $a < |u_R|b$ and its negation. In this last case, and even if it means increasing the value of $a$, once again we will assume that $a > |u^*|b$.

In both cases, (3.9) and (3.10) are still valid, but the definitions of $z^*$ and $\Omega^*$ in (3.11) and (3.13) are modified. More precisely, imposing again the continuity of the Riemann invariants across each wave, easy calculations show that in the case $a < |u_L|b$ and $a < |u_R|b$ formulas in (3.11) are replaced by

$$\left\{ \begin{aligned} z^* &= \frac{|u_R|\tau_R z_R + |u_L|\tau_L z_L}{|u_R|\tau_R + |u_L|\tau_L} - \frac{\Omega_R - \Omega_L}{b|u_R|\tau_R + b|u_L|\tau_L} \\ \Omega^* &= \frac{\Omega_R + \Omega_L}{2} + \frac{b|u_R|\tau_R}{2}(z^* - z_R) - \frac{b|u_L|\tau_L}{2}(z^* - z_L), \end{aligned} \right. \qquad (3.18)$$

while, (3.13) now reads

$$\left\{ \begin{aligned} (z|u|)^* &= |u^*| \frac{\tau_R^* z_R + \tau_L^* z_L}{\tau_R^* + \tau_L^*} - \frac{\Omega_R - \Omega_L}{b(\tau_R^* + \tau_L^*)} \\ \Omega^* &= \frac{\Omega_R + \Omega_L}{2} + \frac{1}{2} (b\tau_R^* ( (z|u|)^* - z_R |u^*| ) - b\tau_L^* ( (z|u|)^* - z_L |u^*| ) ). \end{aligned} \right. \qquad (3.19)$$

At last, instead of using (3.14), here a and b are defined by

$$a = \max(\Delta x, h_L c_L, h_R c_R) \quad \text{and} \quad b = \max\left(\Delta x, \sqrt{g \frac{h_L^2}{u_L^2} \zeta \partial_u (q_b)_L}, \sqrt{g \frac{h_R^2}{u_R^2} \zeta \partial_u (q_b)_R}\right), \qquad (3.20)$$

according to condition (3.17).

## 3.3 Update the bed level in the transport step

For this last strategy we refer to the previous work [13] and we update the bottom height only in the transport step. Thus, the associated acoustic system is now simply given by

$$\left\{ \begin{aligned} \partial_t h + h \partial_x u &= 0 \\ \partial_t (hu) + hu \partial_x u + \partial_x \left(\frac{gh^2}{2}\right) + gh \partial_x z &= 0 \end{aligned} \right. $$
while the transport system has the following form,

$$ egin{cases} 
\partial_t h + u \partial_x h = 0 \\
\partial_t (hu) + u \partial_x (hu) = 0 \\
\partial_t z + \zeta \partial_x q_b = 0.
\end{cases} 
\quad (3.21) $$

Hence, the water height $h$ and flow $hu$ are updated as in the classical shallow-water equations, both in the acoustic and transport steps. In other words, we only need to specify the discretization of the bed level, which does not affect the other two variables. For this reason, this strategy resembles an usual splitting method.

Note that in this context, the relaxation system associated with the acoustic step is nothing but

$$ egin{cases} 
\partial_t \tau - \partial_m u = 0 \\
\partial_t u + \partial_m p = - \frac{g}{\tau} \partial_m z \\
\partial_t \Pi + a^2 \partial_m u = 0
\end{cases} 
\quad (3.22) $$

and the intermediate states associated with the Riemann problems are given by (3.9) and (3.10). At last, let us observe that here the eigenvalues of the relaxation system depend only on the parameter $a$ and not on $b$, with $a$ still defined as in (3.14).

# 4 Numerical method

Before getting into the heart of the matter, we give few details about the time and space discretizations. Given a constant time step $\Delta t$, we define the intermediate times by $t^n = n\Delta t$ for $n \in \mathbb{N}$. Then, the mesh interfaces are $x_{j+1/2} = j\Delta x$ for $j \in \mathbb{Z}$, where $\Delta x$ is the constant space step. Note that for the Lagrangian variable $\xi$, we use the same discretization of the one we introduced for $x$, thus $\Delta x = \Delta \xi$, $x_{j+1/2} = \xi_{j+1/2}$ and $x_j = \xi_j \forall j$, where $x_j$ is the center of the cell $[x_{j-1/2}, x_{j+1/2})$. Hence, given a variable $\varphi$, $\varphi_j^n$ is its piecewise constant approximation at each point $(x_j, t^n)$ with $n \in \mathbb{N}$ and $j \in \mathbb{Z}$. If $\varphi_j^n$ is known, we look for its approximation at the next time level $t^{n+1}$, namely $\varphi_j^{n+1}$. The approximate value obtained at the end of the acoustic step will be denoted by $\varphi_j^{n+1-}$.

Let us note that our numerical schemes are divided into two different steps. First we have the acoustic (Lagrangian) stage, in which we aim at numerically solving one of the systems (3.4), (3.16) or (2.3) depending on the strategy for the bed level equation. Then, we exploit the acoustic (Lagrangian) solution as initial condition for the transport (projection) step, in which we solve either equations (3.2) ($X = h, hu$ and with/without $X = z$) or system (3.21), again depending on the chosen strategy.

## 4.1 Acoustic step

Considering first the acoustic systems (3.3), (3.16) and (2.1), let us first notice that since we have defined approximate Riemann solutions, we can use a classical Godunov-type method, refer to [23, 24, 29] and the references therein. As it is very well-known, it simply consists in averaging on each cell the
juxtaposition of these local approximate Riemann solutions, under a CFL restriction which reads

$$ \Delta t \leq \text{CFL}_l \frac{\Delta x}{\max_j \{ \max(\tau_j^n, \tau_{j+1}^n) \max(a_{j+\frac{1}{2}}, (|u|b)_{j+\frac{1}{2}}) \}}, \quad (4.1) $$

for the acoustic step and

$$ \Delta t \leq \text{CFL}_t \frac{\Delta x}{\max_j \{ u_{j-\frac{1}{2}}^+ - u_{j+\frac{1}{2}}^- \}}, \quad (4.2) $$

for the transport step. CFL$_t$ and CFL$_l$ are respectively the CFL number for the Lagrangian and the transport system, and

$$ u_{j-\frac{1}{2}}^{+} = \max(u_{j-\frac{1}{2}}^{*}, 0) \quad \text{and} \quad u_{j+\frac{1}{2}}^{-} = \min(u_{j+\frac{1}{2}}^{*}, 0). $$

At last, the final time step is taken as the minimum between the two. It is clear that the value of $b$ in (4.1) depends on the numerical treatment of the bottom height $z$ we are employing. Clearly, if $z$ is updated in both steps, we consider formula (3.14) while if the bed elevation is completely taken into account at the acoustic level, $b$ is defined according to (3.20). In the case of the bed level updated only in the transport step, the acoustic time step definition depends only on $a$ and not on $b$.

Regarding the first two equations which are common to (3.3), (3.16) and (2.1), easy calculations not reported here give

$$ \left\{ \begin{aligned} \tau_j^{n+1-} &= \tau_j^n + \frac{\Delta t}{\Delta m_j} (u_{j+\frac{1}{2}}^* - u_{j-\frac{1}{2}}^*) \\ u_j^{n+1-} &= u_j^n - \frac{\Delta t}{\Delta m_j} (\Pi_{j+\frac{1}{2}}^* - \Pi_{j-\frac{1}{2}}^*) - \Delta t \{ \frac{g}{\tau} \partial_m z \}^n_j \end{aligned} \right. \quad (4.3) $$

where $\Delta m_{j+1/2} = (\Delta m_j + \Delta m_{j+1})/2$, $\Delta m_j = \frac{\tau_j^n}{\Delta x}$ and for all $j$

$$ \{\frac{g}{\tau} \partial_m z\}_j^n = \frac{1}{2} \left( \frac{\Delta m_{j+1/2}}{\Delta m_j} \{\frac{g}{\tau} \partial_m z\}_{j+1/2}^n + \frac{\Delta m_{j-1/2}}{\Delta m_j} \{\frac{g}{\tau} \partial_m z\}_{j-1/2}^n \right) \quad \text{with} \quad \{\frac{g}{\tau} \partial_m z\}_{j+1/2}^n = \frac{M_{j+1/2}^n}{\Delta m_{j+1/2}}, $$

the star values $u_{j+\frac{1}{2}}^*$, $\Pi_{j+\frac{1}{2}}^*$ and $M_{j+1/2}$ being locally defined at each interface $x_{j+1/2}$ and for all $j$ thanks to (3.9) and (3.10). Here of course, the subscripts $L$ and $R$ stand for $j$ and $j+1$ respectively.

In Lagrangian coordinates, (4.3) reveals to be strictly equivalent to

$$ \left\{ \begin{aligned} L_j^{n+1-} h_j^{n+1-} &= L_j^n h_j^n \\ L_j^{n+1-} (h u)_j^{n+1-} &= L_j^n (h u)_j^n - \frac{\Delta t}{\Delta x} (\Pi_{j+\frac{1}{2}}^* - \Pi_{j-\frac{1}{2}}^*) - \Delta t \{ g h \partial_x z \}_j^n \end{aligned} \right. \quad (4.4) $$

where

$$ L_j^{n+1-} = L_j^n + \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* - u_{j-\frac{1}{2}}^*) \quad \text{with} \quad L_j^n = 1 $$

and

$$ s = -gh\partial_x z, \quad s_j^n = \frac{1}{2} (s_{j+1/2}^n + s_{j-1/2}^n), \quad s_{j+1/2}^n = -\frac{M_{j+1/2}^n}{\Delta x} \quad \forall j. $$
Let us now briefly give the update formulas for the topography.

*Bed level in both steps.* Considering system (3.4), we state

$$L_j^{n+1-} z_j^{n+1-} = L_j^n z_j^n - \frac{\Delta t}{\Delta x} ((\Omega - zu)_{j+\frac{1}{2}}^* - (\Omega - zu)_{j-\frac{1}{2}}^*)$$

which turns out to be equivalent to

$$z_j^{n+1-} (1 + \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* - u_{j-\frac{1}{2}}^*)) = z_j^n + \frac{\Delta t}{\Delta x} ((zu)_{j+\frac{1}{2}}^* - (zu)_{j-\frac{1}{2}}^*) - \frac{\Delta t}{\Delta x} (\Omega_{j+\frac{1}{2}}^* - \Omega_{j-\frac{1}{2}}^*)$$

where the star values are easily defined from (3.11) and (3.13).

*Bed level in the acoustic step.* In this case, from system (3.15) we find

$$z_j^{n+1-} = z_j^n - \frac{\Delta t}{\Delta x} (\Omega_{j+\frac{1}{2}}^* - \Omega_{j-\frac{1}{2}}^*), \quad (4.5)$$

where $\Omega_{j\pm\frac{1}{2}}^*$ are defined from (3.18) and (3.19).

*Bed level in the transport step.* Finally, it is clear that here we simply have $z_j^{n+1-} = z_j^n$ as the bed level is completely taken into account in the transport step.

## 4.2 Transport step

Referring to [37], we can approximate $\partial_t X + u\partial_x X = 0$ by

$$X_j^{n+1} = (LX)_j^{n+1-} - \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* (LX)_{j+\frac{1}{2}}^{n+1-} - u_{j-\frac{1}{2}}^* (LX)_{j-\frac{1}{2}}^{n+1-}), \quad (4.6)$$

where

$$(LX)_{j+\frac{1}{2}}^{n+1} = \begin{cases} (LX)_j^{n+1-} & \text{if } u_{j+\frac{1}{2}}^* \ge 0 \\ (LX)_{j+1}^{n+1-} & \text{if } u_{j+\frac{1}{2}}^* < 0 \end{cases}$$

with either $X = h$ or $X = hu$. Notice that this formula can also be explained in an alternative way, using the Lagrangian coordinates. Indeed, defining $\hat{\xi}_{j+\frac{1}{2}}(t)$ such that for all $j$

$$x(\hat{\xi}_{j+\frac{1}{2}}(T), T) = x_{j+\frac{1}{2}}, \quad \text{with } T \ge 0,$$

it is enough to recall that

$$X_j(t) = \frac{1}{\Delta x} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} X(x, t) dx = \frac{1}{\Delta x} \int_{x(\hat{\xi}_{j-\frac{1}{2}}, t)}^{x(\hat{\xi}_{j+\frac{1}{2}}, t)} X(x, t) dx = \frac{1}{\Delta x} \int_{\hat{\xi}_{j-\frac{1}{2}}}^{\hat{\xi}_{j+\frac{1}{2}}} L(\xi, t) \bar{X}(\xi, t) d\xi$$
and then split the last integral into three parts to define $X_j^{n+1}$, namely

$$X_j^{n+1} = \frac{1}{\Delta x} \int_{\hat{\xi}_{j-\frac{1}{2}}}^{\xi_{j-\frac{1}{2}}} L(\xi, t^{n+1-}) \bar{X}(\xi, t^{n+1-}) d\xi + \\ + \frac{1}{\Delta x} \int_{\xi_{j-\frac{1}{2}}}^{\xi_{j+\frac{1}{2}}} L(\xi, t^{n+1-}) \bar{X}(\xi, t^{n+1-}) d\xi + \frac{1}{\Delta x} \int_{\xi_{j+\frac{1}{2}}}^{\hat{\xi}_{j+\frac{1}{2}}} L(\xi, t^{n+1-}) \bar{X}(\xi, t^{n+1-}) d\xi \quad (4.7)$$

where we approximate $x_{j+\frac{1}{2}}$ at first-order, $x_{j+\frac{1}{2}} = x(\hat{\xi}_{j+\frac{1}{2}}(T), T) \approx x(\hat{\xi}_{j+\frac{1}{2}}(T), 0) + T \partial_t x(\hat{\xi}_{j+\frac{1}{2}}(T), 0) \approx \hat{\xi}_{j+\frac{1}{2}} + T u_{j+\frac{1}{2}}^*$, for a fixed time $T \ge 0$. Then, it is clear that formula (4.6) can be seen as a first-order approximation of the integrals in (4.7). Note that the second-order scheme will be obtained by approximating the three integrals in (4.7) at second-order of accuracy in space. Since this procedure has been explained in details in [37], here we do not provide further information.

As before, let us now give the update formulas for the topography.

**Bed level in both steps.** In this case, similarly to the formulas for $h, hu$, we have for all $j$

$$z_j^{n+1} = (Lz)_j^{n+1-} - \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* z_{j+\frac{1}{2}}^{n+1-} - u_{j-\frac{1}{2}}^* z_{j-\frac{1}{2}}^{n+1-})$$

with the only difference that we use the values of $z_{j+\frac{1}{2}}^{n+1-}$ and not of its Lagrangian counterpart.

**Bed level in the acoustic step.** Since in this numerical method the bottom height is completely considered in the acoustic step, here we simply have

$$z_j^{n+1} = z_j^{n+1-}.$$

**Bed level in the transport step.** In this case, it is a matter of discretizing the full Exner equation, which can be done following at least two options. Referring to the previous work [13], a first option called *decoupled* consists in simply updating the topography as in the following

$$z_j^{n+1} = z_j^n - \zeta \frac{\Delta t}{\Delta x} \left( u_{j+\frac{1}{2}}^* \left( \frac{q_b(u)}{u} \right)_{j+\frac{1}{2}}^n - u_{j-\frac{1}{2}}^* \left( \frac{q_b(u)}{u} \right)_{j-\frac{1}{2}}^n \right), \quad (4.8)$$

with

$$\left(\frac{q_b(u)}{u}\right)_{j+1/2}^n = \begin{cases} \left(\frac{q_b(u)}{u}\right)\left(h_{j+1}^n, (hu)_j^n\right) & \text{if } u_{j+1/2}^* \le 0 \\ \left(\frac{q_b(u)}{u}\right)\left(h_j^n, (hu)_j^n\right) & \text{if } u_{j+1/2}^* > 0. \end{cases} \quad (4.9)$$

Note that $(\frac{q_b(u)}{u})_{j+1/2}^n$ can be defined with no ambiguity since we either suppose the solid transport discharge to be given by the Grass formula (1.2) with $m_g = 3$ (and therefore $u$ simplifies) or by the MPM formula (1.3) where $q_b$ is considered null for $u = 0$. It is clear that this updating formula depends only on the solution at the initial time, and not on the solution obtained at the end of the acoustic step, which justifies the term *decoupled*.
As a second option called weakly coupled, one can also exploit the values of Lu obtained at the end of the acoustic step and thus set

$$z_j^{n+1} = z_j^n - \zeta \frac{\Delta t}{\Delta x} \left( u_{j+\frac{1}{2}}^* \left( \frac{q_b(Lu)}{Lu} \right)_{j+\frac{1}{2}}^{n+1-} - u_{j-\frac{1}{2}}^* \left( \frac{q_b(Lu)}{Lu} \right)_{j-\frac{1}{2}}^{n+1-} \right). \quad (4.10)$$

and $(\frac{q_b(Lu)}{Lu})_{j+\frac{1}{2}}^{n+1-}$ defined as the corresponding value in the $(x_j, t^{n+1-})$ or $(x_{j+1}, t^{n+1-})$ point respectively if $u_{j+\frac{1}{2}}^* > 0$ or $u_{j+\frac{1}{2}}^* \le 0$ and according to (4.9). Since a definition of $(Lu)^{n+1-}$ is needed, we observe that an evolution equation for Lu in the Lagrangian step reads

$$\partial_t (Lu) - \partial_\xi \frac{u^2}{2} = -g \partial_\xi (h + z),$$

and can be discretized as

$$(Lu)_j^{n+1-} = (Lu)_j^n + \frac{\Delta t}{2\Delta x} ((u_{j+\frac{1}{2}}^*)^2 - (u_{j-\frac{1}{2}}^*)^2) + \Delta t \frac{\hat{s}_{j+\frac{1}{2}} + \hat{s}_{j-\frac{1}{2}}}{2} \quad (4.11)$$

where $\hat{s}_{j+\frac{1}{2}} = -g((h+z)_{j+1} - (h+z)_j)/\Delta x$, which concludes the definition of the schemes.

## 4.3 Friction term approximation

As far as the friction term approximation is concerned, here we refer to the work of Audusse et al. [3]. Hence, we briefly recall the discretization they used and refer to their paper for more details. In particular, we exploit an implicit splitting strategy. Once the solution $hu^{n+1}$ from the projection step has been obtained, we state

$$\bar{h}u^{n+1} = hu^{n+1} - g\Delta t \frac{\mu_j^2 |\bar{h}u^{n+1}| \bar{h}u^{n+1}}{h^{n+1} R_h^{4/3}} \quad (4.12)$$

which gives us the flux $hu$ at the new time level with the friction contributions included. Note that, imposing $\bar{h}u^{n+1}$ and $hu^{n+1}$ to have the same sign, it is possible to obtain the explicit solution of equation (4.12) so that the computational cost is not high. Moreover, as explained in [3], this discretization will also allow us to preserve the "constant bed slope" equilibrium (1.5).

# 5 Increasing the order of accuracy

So far we presented three different first-order numerical schemes for the shallow water Exner system. The discretizations for the water height and flow appeared to be the same in all the methods, only the manner in which we updated the bottom height changed. We now aim to develop second-order numerical methods, which could lead to the construction of even higher order schemes. As we will see, we do not achieve the second-order of accuracy for all the three schemes in the same way. This is mainly due to the fact that we also ask for the well-balanced property and it cannot be obtained in the same manner as it depends on the underlying first-order scheme. However, the heart of the method is
the same, namely we exploit polynomial reconstructions [41] and Runge-Kutta TVD approach [26] in order to reach the second order of accuracy in space and time respectively. Indeed, what will change is mainly the way we define the slopes for the reconstruction polynomials in space as the Runge-Kutta procedure does not affect the ability of the scheme of preserving the well-balanced property. We specify that the Runge-Kutta scheme is used at second order and applied to the overall scheme, namely the acoustic and transport step together.

## 5.1 Update the bed level in both steps

Let us start with the numerical scheme which entails a splitting of the Exner equation in both steps. Regarding the acoustic step, we first proceed in a very classical way by making use of conservative and first-order polynomial reconstructions of the form

$$ \mathbf{P}_j(x) = \mathbf{Q}_j^n + \Delta_j(x - x_j), \qquad (5.1) $$

where $\mathbf{Q} = (h, hu, z)^t$ is the vector of unknowns and $\Delta_j = (\Delta_j(h), \Delta_j(hu), \Delta_j(z))^t$ denotes the corresponding slopes. Motivated by the well-balanced property, we compute the slopes using standard ENO or MINMOD limiters applied to the free surface $H = h+z, hu$ and $z$. Then, we simply set $\Delta_j(h) = \Delta_j(H) - \Delta_j(z)$. Indeed, notice that $H$ is constant under the "lake at rest" condition and therefore the slopes $\Delta_j(H)$ automatically reduce to zero in this case. Then, it clearly follows $\Delta_j(h) = -\Delta_j(z)$, which is necessary for the well-balanced property.

Finally, the Lagrangian step reads

$$ \left\{ \begin{aligned} L_j^{n+1-} h_j^{n+1-} &= L_j^n h_j^n \\ L_j^{n+1-} (hu)_j^{n+1-} &= L_j^n (hu)_j^n - \frac{\Delta t}{\Delta x} (\Pi_{j+\frac{1}{2}}^* - \Pi_{j-\frac{1}{2}}^*) + \Delta t (s_j^n + s_{C,j}^n) \\ L_j^{n+1-} z_j^{n+1-} &= L_j^n z_j^n - \frac{\Delta t}{\Delta x} ((\Omega - zu)^*_{j+\frac{1}{2}} - (\Omega - zu)^*_{j-\frac{1}{2}}) \end{aligned} \right. \qquad (5.2) $$

where we compute the interface values using the first-order formulas but applied to the left and right traces of the reconstruction polynomials, namely

$$ \mathbf{Q}_{j+\frac{1}{2}L} = \mathbf{P}_j(x_{j+\frac{1}{2}}) \quad \text{and} \quad \mathbf{Q}_{j+\frac{1}{2}R} = \mathbf{P}_{j+1}(x_{j+\frac{1}{2}}), \qquad (5.3) $$

instead of $\mathbf{Q}_j$ and $\mathbf{Q}_{j+1}$. Notice that an additional term $s_{C,j}^n$ defined by

$$ s_{C,j}^n = -g \frac{h_{j-\frac{1}{2}R} + h_{j+\frac{1}{2}L}}{2} \frac{z_{j+\frac{1}{2}L} - z_{j-\frac{1}{2}R}}{\Delta x} $$

and representing the in-cell second-order contribution of the source term is introduced to make the scheme well-balanced. Of course, we note that $s_{C,j}^n = 0$ for all $j$ if the slopes are null.

Regarding the transport step, we exploit again polynomial reconstructions but now we first reconstruct the Lagrangian variables $(Lh)_j^{n+1-}$ and $(Lhu)_j^{n+1-}$ obtained at the end of the Lagrangian step. Then, the updating formula for $h$ and $hu$ are given by a second-order approximation of the three integrals that appear in (4.7). This is achieved by using a classical mid-point rule. As far as the topography
is concerned, the procedure is similar, the only difference is that we reconstruct its values $z_j^{n+1-}$ instead of $(Lz)_j^{n+1-}$. Hence, the reconstructed polynomial is given by

$$P_j(z, \xi) = z_j^{n+1-} + \Delta_j(\xi - \xi_j),$$

and the second-order updating formula for $z$ simply reads

$$z_j^{n+1} = (Lz)_j^{n+1-} - \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* P_{j+\frac{1}{2}} \left( z, \frac{\xi_{j+\frac{1}{2}} + \hat{\xi}_{j+\frac{1}{2}}}{2} \right) - u_{j-\frac{1}{2}}^* P_{j-\frac{1}{2}} \left( z, \frac{\xi_{j-\frac{1}{2}} + \hat{\xi}_{j-\frac{1}{2}}}{2} \right)),$$

with

$$P_{j-\frac{1}{2}}(z, \xi) = \begin{cases} P_{j-1}(z, \xi) & \text{if } \xi_{j-\frac{1}{2}} > \hat{\xi}_{j-\frac{1}{2}} \\ P_j(z, \xi) & \text{if } \xi_{j-\frac{1}{2}} \le \hat{\xi}_{j-\frac{1}{2}}. \end{cases}$$

## 5.2 Update the bed level in the acoustic step

Let us proceed with the numerical scheme that takes into account the bed level $z$ only at the acoustic level. Here, we follow exactly the same above-described procedure to obtain a second-order discretization for the variables $h$ and $hu$. The only difference is related to the bottom height approximation for which we use formula (4.5) and either (3.18) or (3.19) for the star values. We underline that, in the acoustic step, the same reconstruction procedure is considered for the variables $h, hu$ and $z$ whereas, in the transport step, nothing has to be done for $z$.

## 5.3 Update the bed level in the transport step

For this last strategy we refer again to [13] and we give few details. The Lagrangian step simply reads

$$\begin{cases} L_j^{n+1-} h_j^{n+1-} = L_j^n h_j^n \\ L_j^{n+1-} (hu)_j^{n+1-} = L_j^n (hu)_j^n - \frac{\Delta t}{\Delta x} (\Pi_{j+\frac{1}{2}}^* - \Pi_{j-\frac{1}{2}}^*) + \Delta t s_j^n, \end{cases}$$

where again the star values are computed using left and right traces of piecewise linear approximations of the solution at each interface, instead of piecewise constant approximation like in the first-order method. Here, and in order to maintain the well-balanced property, we exploit the so-called fluctuations [37, 12] to define the slopes of the reconstructed polynomials. In particular, notice that here the bed level is kept constant in the Lagrangian step and therefore no reconstruction is applied to this variable (the corresponding slopes are considered to be zero). As far as the transport step is concerned, the water height $h$ and the discharge $hu$ are updated like in the two other methods and thus we do not further discuss it. Regarding the topography, we consider (4.10) and, in order to make it second order in space, in place of the piecewise constant approximation of $Lu$ given by (4.11) we use its polynomial reconstructed values.
# 6 Overall discretization and well-balanced property

This section is devoted to the illustration of the first-order overall discretizations of the previous schemes and their well-balanced property. Here we also take advantage of this section to observe that the three first-order methods we presented are also positivity-preserving under suitable conditions.

**Remark 1.** Is it possible to prove that the three first-order methods are also positivity-preserving, meaning that they are able to preserve the strict positivity of the water height under the associated CFL conditions $4.1$ and $4.2$ with $CFL_t \le \frac{1}{2}$ and $CFL_t \le 1$. See also [11] for more details.

## 6.1 Bed level in both steps

Considering the numerical treatment which update $z$ in both the acoustic and transport step (sections 3.1, 4 and 5.1), it is easy to see that the whole scheme at first-order takes the following final form,

$$ egin{cases} h_j^{n+1} = h_j^n - \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^*(Lh)_{j+\frac{1}{2}}^{n+1-} - u_{j-\frac{1}{2}}^*(Lh)_{j-\frac{1}{2}}^{n+1-}) \\ (hu)_j^{n+1} = (hu)_j^n - \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^*(Lhu)_{j+\frac{1}{2}}^{n+1-} + \Pi_{j+\frac{1}{2}}^* - (u_{j-\frac{1}{2}}^*(Lhu)_{j-\frac{1}{2}}^{n+1-} + \Pi_{j-\frac{1}{2}}^*)) + \Delta t s_j^n, \\ z_j^{n+1} = z_j^n - \frac{\Delta t}{\Delta x} ((\Omega - zu)^*_{j+\frac{1}{2}} - (\Omega - zu)^*_{j-\frac{1}{2}}) - \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* z_{j+\frac{1}{2}}^{n+1-} - u_{j-\frac{1}{2}}^* z_{j-\frac{1}{2}}^{n+1-}) \end{cases} \quad (6.1) $$

with

$$ X_{j+\frac{1}{2}}^{n+1-} = \begin{cases} X_j^{n+1-} & \text{if } u_{j+\frac{1}{2}}^{*,n+1-} > 0 \\ X_{j+1}^{n+1-} & \text{if } u_{j+\frac{1}{2}}^{*,n+1-} \le 0 \end{cases} \quad (6.2) $$

and $X = Lh, Lhu, z$. Note that the evolution equation for the topography $z$ can also be reformulated as

$$ z_j^{n+1} = z_j^n - \frac{\Delta t}{\Delta x} (\Omega_{j+\frac{1}{2}}^* - \Omega_{j-\frac{1}{2}}^*) + \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* (z_{j+\frac{1}{2}}^* - z_{j+\frac{1}{2}}^{n+1-}) - u_{j-\frac{1}{2}}^* (z_{j-\frac{1}{2}}^* - z_{j-\frac{1}{2}}^{n+1-})). $$

It is clear that without the source term present in the evolution equation for $hu$, the whole numerical scheme would be conservative.

Let us now prove that both the first and second-order schemes are well-balanced.

**Theorem 1.** The numerical method with updating formula (6.1) and star values given in section 3.1.1 is well-balanced under the "lake at rest" condition (1.4).

*Proof.* Referring to [16], it can be easily seen that $u_{j+1/2}^* = 0$, $h_j^{n+1-} = h_j^n$ and $(hu)_j^{n+1-} = (hu)_j^n \forall j$. Let us now consider the topography $z$. Since we are under the hypothesis that $u_j^n = u_{j+1/2}^* = 0 \forall j$, it appears clear that we are in the case $|u_{j+1/2}^*|b_{j+1/2} < a_{j+1/2}$. Consequently we should use formulae (3.13) in order to update $z_j^{n+1-}$, which leads us to $(|u|z)_{j+1/2}^* = 0$ and $\Omega_{j+1/2}^* = 0 \forall j$. Thus the scheme is well-balanced for $z$ as well. □

**Theorem 2.** The second-order numerical method which update the bed level $z$ in both the acoustic and transport step (formulas in section 5.1) preserves the "lake at rest" stationary solution (1.4).
Proof. First of all, we observe that $u_{j+\frac{1}{2}}^* = 0 \forall j$ when the variables at time $t^n$ satisfy the "lake at rest" solution. Indeed, we have

$$u_{j+\frac{1}{2}}^* = \frac{1}{2}(u_{j+\frac{1}{2}L} + u_{j+\frac{1}{2}R}) - \frac{1}{2a_{j+\frac{1}{2}}}(\Pi_{j+\frac{1}{2}R} - \Pi_{j+\frac{1}{2}L}) - \frac{1}{2a_{j+\frac{1}{2}}}(\frac{g}{2}(h_{j+\frac{1}{2}L} + h_{j+\frac{1}{2}R})(z_{j+\frac{1}{2}R} - z_{j+\frac{1}{2}L}))$$

but $\Delta_j^t(u) = 0 \forall j$ and thus $u_{j+\frac{1}{2}L} = u_{j+\frac{1}{2}R} = u_j = u_{j+1} = 0$. Thus, we only have to prove that

$$\Pi_{j+\frac{1}{2}R} - \Pi_{j+\frac{1}{2}L} = -\frac{g}{2}(h_{j+\frac{1}{2}L} + h_{j+\frac{1}{2}R})(z_{j+\frac{1}{2}R} - z_{j+\frac{1}{2}L}) \quad (6.3)$$

but

$$\Pi_{j+\frac{1}{2}R} - \Pi_{j+\frac{1}{2}L} = \frac{g}{2}(h_{j+\frac{1}{2}R}^2 - h_{j+\frac{1}{2}L}^2) = \frac{g}{2}(h_{j+\frac{1}{2}L} + h_{j+\frac{1}{2}R})(h_{j+\frac{1}{2}R} - h_{j+\frac{1}{2}L})$$

and finally (6.3) is equivalent to

$$h_{j+\frac{1}{2}R} - h_{j+\frac{1}{2}L} = -(z_{j+\frac{1}{2}R} - z_{j+\frac{1}{2}L}).$$

Exploiting definitions (5.1)-(5.3), we write

$$\begin{aligned} h_{j+\frac{1}{2}R} - h_{j+\frac{1}{2}L} &= h_{j+1} - \frac{\Delta_{j+1}^t(h)\Delta x}{2} - (h_j + \frac{\Delta_j^t(h)\Delta x}{2}) = \\ &\quad -z_{j+1} + z_j + \frac{\Delta_{j+1}^t(z)\Delta x}{2} + \frac{\Delta_j^t(z)\Delta x}{2} = -z_{j+\frac{1}{2}R} + z_{j+\frac{1}{2}L}. \end{aligned}$$

Thus, we proved that $u_{j+\frac{1}{2}}^* = 0 \forall j$. Next, if the "lake at rest" condition holds true, from the second equation of system (5.2) we observe

$$-(\Pi_{j+\frac{1}{2}}^* - \Pi_{j-\frac{1}{2}}^*) + \Delta x(s_j^n + s_{C,j}^n) = 0$$

and thus $L_j^{n+1-}(hu)_j^{n+1-} = L_j^n(hu)_j^n = 0$. Let us give the details. Since we know $u_{j\pm\frac{1}{2}L} = u_{j\pm\frac{1}{2}R} = 0 \forall j$, we can write

$$\begin{aligned} \Pi_{j+\frac{1}{2}}^* - \Pi_{j-\frac{1}{2}}^* &= \frac{1}{2}(\Pi_{j+\frac{1}{2}R} + \Pi_{j+\frac{1}{2}L}) - \frac{1}{2}(\Pi_{j-\frac{1}{2}R} + \Pi_{j-\frac{1}{2}L}) = \\ &= \frac{1}{2}((\Pi_{j+\frac{1}{2}R} - \Pi_{j+\frac{1}{2}L}) + (\Pi_{j-\frac{1}{2}R} - \Pi_{j-\frac{1}{2}L})) + \Pi_{j+\frac{1}{2}L} - \Pi_{j-\frac{1}{2}R} \\ \stackrel{(6.3)}{=} \Delta x s_j^n + \Pi_{j+\frac{1}{2}L} - \Pi_{j-\frac{1}{2}R}. \end{aligned}$$

With similar computations, it is straightforward to see that $\Pi_{j+\frac{1}{2}L} - \Pi_{j-\frac{1}{2}R} = \Delta x s_{C,j}^n$. Finally, since $u_{j\pm\frac{1}{2}}^* = 0 \forall j$, we also know that $\Omega_{j\pm\frac{1}{2}}^* = 0 \forall j$ and thus $z_j^{n+1-} = z_j^n$.

Regarding the transport step, it does not modify the unknowns' values as $u_{j\pm\frac{1}{2}}^* = 0 \forall j$. Similarly, since we already know that the Runge-Kutta procedure does not prevent the scheme to be well-balanced, the property is proved. $\square$
**Theorem 3.** The numerical method with updating formula (6.1) and star values given in section 3.1.1 is well-balanced under the "constant bed slope" equilibrium (1.5). This statement remains true even if we consider its second-order accurate version presented in section 5.1.

*Proof.* Assuming to be under the "constant bed slope" equilibrium (1.5), namely to have constant velocity $u_i^n = u_{i+1}^n \forall i$, constant water height $h_i^n = h_{i+1}^n \forall i$ and constant bed slope $z_{i+1}^n - z_i^n = -\Delta x S_{f,i} = K \forall i$ with constant $K$, we want to prove that such a steady state is preserved, namely that $h_i^{n+1} = h_i^n$, $hu_i^{n+1} = hu_i^n$ and $z_i^{n+1} = z_i^n \forall i$. Let us start seeing what are the star values in this case. We consider only the case $a < |u_L|b, a < |u_R|b$ as for the other one the procedure is analogous. Thus, it easy to see that we obtain $\Pi_{i+1/2}^* = \Pi_i^n = \Pi^* = \text{constant}$ and $u_{i+1/2}^* = u_i^n - \frac{1}{2a}M_{i+1/2} = u^* = \text{constant}$ as also the parameter $a$ and $b$ are clearly constant in all the domain. Moreover, we also have $z_{i+1/2}^* = \frac{z_{i+1}^n + z_i^n}{2} + \frac{u}{2|u|b\tau}(z_{i+1}^n - z_i^n)$ and $\Omega_{i+1/2}^* = \Omega^* = \text{constant}$. Then, it is easy to see that we obtain $h_i^{n+1} = Lh_i^{n+1-} = h_i^n$, while for the discharge we get

$$hu_i^{n+1} = hu_i^n + \Delta tgh_i^n S_{f,i}$$

where the friction has not yet been taken into account and which compensates the presence of the term $\Delta tgh_i^n \Delta x S_f$, leading to $\bar{h}u_i^{n+1} = hu_i^n$, see also [3]. Finally, as far as the topography is concerned, we have

$$ \begin{aligned} z_i^{n+1} &= z_i^n - \frac{\Delta t}{\Delta x}(\Omega_{i+1/2}^* - \Omega_{i-1/2}^*) + \frac{\Delta t}{\Delta x}(u_{i+1/2}^*(z_{i+1/2}^* - z_{i+1/2}^{n+1-}) - u_{i-1/2}^*(z_{i-1/2}^* - z_{i-1/2}^{n+1-})) \\ &= z_i^n + u_i^* \frac{\Delta t}{\Delta x}((z_{i+1/2}^* - z_{i+1/2}^{n+1-}) - (z_{i-1/2}^* - z_{i-1/2}^{n+1-})). \end{aligned} $$

Then, since $u^*$ is constant, it is clear that we have $z_{i+1/2}^{n+1-} - z_{i-1/2}^{n+1-} = z_{i+1}^{n+1-} - z_i^{n+1-} = z_i^{n+1-} - z_{i-1}^{n+1-} = z_i^n - z_i^n = z_i^n - z_{i-1}^n = K$ and thus

$$z_i^{n+1} = z_i^n + u_i^* \frac{\Delta t}{\Delta x} \left( \frac{z_{i+1}^n + z_i^n}{2} - \frac{z_i^n + z_{i-1}^n}{2} - K \right) = z_i^n + u_i^* \frac{\Delta t}{\Delta x} \left( \frac{z_{i+1}^n - z_i^n}{2} + \frac{z_i^n - z_{i-1}^n}{2} - K \right) = z_i^n,$$

which concludes the proof for the first-order method. Concerning the second-order extension, it is enough to check that the slopes of the polynomial reconstruction are either null or constant and similar computations lead to the same result; thus we do not insert further details here. $\square$

## 6.2 Bed level in the acoustic step

The first-order numerical method, which takes into account the bed level $z$ only in the acoustic step (sections 3.2, 4 and 5.2), has an overall form similar to the one (6.1) seen in the previous section, the only difference related to the update equation for $z$. In practice, we write

$$ \left\{ \begin{aligned} h_j^{n+1} &= h_j^n - \frac{\Delta t}{\Delta x} (u_{j+1/2}^* L h_{j+1/2}^{n+1-} - u_{j-1/2}^* L h_{j-1/2}^{n+1-}) \\ hu_j^{n+1} &= hu_j^n - \frac{\Delta t}{\Delta x} (u_{j+1/2}^* L h u_{j+1/2}^{n+1-} + \Pi_{j+1/2}^* - (u_{j-1/2}^* L h u_{j-1/2}^{n+1-} + \Pi_{j-1/2}^*)) + \Delta t s_j^n, \end{aligned} \right. \quad (6.4) $$

with $X_{j+1/2}^{n+1-}$ given by formula (6.2) and $X = Lh, Lhu$.
**Theorem 4.** The numerical method with updating formula (6.4) and star values given in section 3.2.2 is well-balanced under the "lake at rest" condition (1.4).

*Proof.* Proof analogous to the one of theorem 1.

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

**Theorem 5.** The second-order numerical method which update the bed level $z$ only in the acoustic step (section 5.2) preserves the "lake at rest" stationary solution (1.4).

*Proof.* Proof analogous to the one of theorem 2.

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

**Theorem 6.** The numerical method with updating formula (6.4) and star values given in section 3.2.2 is well-balanced under the "constant bed slope" equilibrium (1.5). This statement remains true even if we consider its second-order accurate version presented in section 5.2.

*Proof.* Proof analogous to the one of theorem 3.

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

## 6.3 Bed level in the transport step

Finally, let us briefly comment on the numerical scheme which resembles to the usual splitting method (sections 3.3, 4 and 5.3). The bed level is updated only in the transport step and, as such, the first-order overall discretization reads

$$ egin{cases} h_j^{n+1} = h_j^n - \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* L h_{j+\frac{1}{2}}^{n+1-} - u_{j-\frac{1}{2}}^* L h_{j-\frac{1}{2}}^{n+1-}) \\ hu_j^{n+1} = hu_j^n - \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* L h_{j+\frac{1}{2}}^{n+1-} + \Pi_{j+\frac{1}{2}}^* - (u_{j-\frac{1}{2}}^* L h_{j-\frac{1}{2}}^{n+1-} + \Pi_{j-\frac{1}{2}}^*)) + \Delta t s_j^n \\ z_j^{n+1} = z_j^n - \zeta \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* \left( \frac{q_b(Lu)}{Lu} \right)^{n+1-}_{j+\frac{1}{2}} - u_{j-\frac{1}{2}}^* \left( \frac{q_b(Lu)}{Lu} \right)^{n+1-}_{j-\frac{1}{2}}), \end{cases} \quad (6.5) $$

see formulas (4.4), (4.6) and (4.10). Referring to [13], we state that both the first and second-order accurate versions of this method are well-balanced.

**Theorem 7.** The first-order numerical scheme with updating formula (6.5) preserves the "lake at rest" stationary solution (1.4).

*Proof.* Refer to the previous work [13].

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

**Theorem 8.** The second-order numerical method which update the bed level $z$ only in the transport step (section 5.3) preserves the "lake at rest" stationary solution (1.4).

*Proof.* Refer to the previous work [13].

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

**Theorem 9.** The first-order numerical scheme with updating formula (6.5) preserves the "constant bed slope" equilibrium (1.5). This statement remains true even if we consider its second-order accurate version presented in section 5.3.

*Proof.* Proof analogous to the one of theorem 3.

<figure><img src="image_7.png" /></figure>
# 7 Numerical evidences

Here we test the numerical schemes we presented so far; for the sake of conciseness we distinguish them by calling them as in the following:

* "AcTrZ": scheme with bed level $z$ updated in both the ACoustic and the TRansport step (see sections 3.1, 4 and 5.1);

* "AcZ": scheme with bed level $z$ updated in the ACoustic step (see sections 3.2, 4 and 5.2);

* "TrZ": scheme with bed level $z$ updated in the TRansport step (see sections 3.3, 4 and 5.3), using in particular discretization (4.10).

If not otherwise specified, we take $\zeta = 1$, $q_b = A_g u^3$ with $A_g = 0.005$ for the Exner equation, $\mu_f = 0$ and transmissive boundary conditions. For the CFL number we use $CFL_t = 0.45$ and $CFL_l = 0.25$ for the first and second order schemes respectively, while $CFL_t = 1$. With the exception of the accuracy test case, when the reference solution is inserted, it is computed exploiting the second-order AcTrZ scheme with $M = 1000$ cells, where $\Delta x = \frac{L}{M}$ with $L$ the length of the channel.

## 7.1 Test of order of accuracy

Here we test the order of accuracy of the numerical schemes described previously. Let us consider a channel of length $L = 20m$, $A_g = 0.3$, $m = 3$. The initial condition is given by null velocity and

$$ \begin{cases} z_{IC} = 0.1 - 0.01e^{-(x-10)^2} \\ h_{IC} = 2 - 0.1e^{-(x-10)^2}. \end{cases} $$

We refer to paper [7] for this test case. The reference solution is computed using $M = 2048$ cells and TrZ-second-order "decoupled" method (formula (4.8)). In table 1 we insert the errors and the EOA in norm $L^1$ for all the three schemes. We can see that the second-order of accuracy is reached in each case.

## 7.2 Riemann problem: dam break on movable bottom

For this Riemann problem we refer to [1]. The length of the channel is $L = 10m$ and the dam is placed in the middle. The ending time is $t_{end} = 1s$. The initial condition is given by null velocity, flat topography and water height $h_L = 2m$ if $x < L/2$, $h_R = 0.125m$ if $x > L/2$. First of all in figure 1 we compare the first and second order method, in particular we use LP-AcZ but similar results can be found with the other two schemes. Clearly, the second-order version of the method gives more accurate results for the same value of the mesh. On the right we used $M = 2000$ cells to show that the first-order solution converges to the second-order one in general, even if we can observe a small difference of order $10^{-2}$ in the shock position. However, this is not surprising but natural, as the shallow water Exner system is not conservative. Then, in figure 2 we insert the bed level solution of all the three second-order numerical schemes in order to be able to sum up the merits and flaws of each of them. As expected, the more diffusive numerical scheme is the one in which the bed level is updated in both steps,
<table><caption>Table 1: Errors and empirical convergence rates for norm L¹. Mesh of size M = (64, 128, 256, 512, 1024), CFL = 0.25.</caption><thead><tr><th>Method</th><th>Mesh M</th><th>Variable</th><th>err L<sup>1</sup></th><th>O(L<sup>1</sup>)</th><th>Variable</th><th>err L<sup>1</sup></th><th>O(L<sup>1</sup>)</th><th>Variable</th><th>err L<sup>1</sup></th><th>O(L<sup>1</sup>)</th></tr></thead><tbody><tr><td rowspan="5">AcTrZ-scheme</td><td>64</td><td rowspan="5">h</td><td>0.0269</td><td>—</td><td rowspan="5">hu</td><td>0.1174</td><td>—</td><td rowspan="5">z</td><td>0.1372 × 10<sup>-3</sup></td><td>—</td></tr><tr><td>128</td><td>0.0083</td><td>1.6927</td><td>0.0354</td><td>1.7294</td><td>0.0639 × 10<sup>-3</sup></td><td>1.1030</td></tr><tr><td>256</td><td>0.0027</td><td>1.6190</td><td>0.0115</td><td>1.6236</td><td>0.0218 × 10<sup>-3</sup></td><td>1.5529</td></tr><tr><td>512</td><td>0.0007</td><td>1.8757</td><td>0.0031</td><td>1.8773</td><td>0.0061 × 10<sup>-3</sup></td><td>1.8456</td></tr><tr><td>1024</td><td>0.0002</td><td>1.9786</td><td>0.0008</td><td>1.9819</td><td>0.0016 × 10<sup>-3</sup></td><td>1.9532</td></tr><tr><td rowspan="5">AcZ-scheme</td><td>64</td><td rowspan="5">h</td><td>0.0269</td><td>—</td><td rowspan="5">hu</td><td>0.1174</td><td>—</td><td rowspan="5">z</td><td>0.1370 × 10<sup>-3</sup></td><td>—</td></tr><tr><td>128</td><td>0.0083</td><td>1.6931</td><td>0.0354</td><td>1.7301</td><td>0.0592 × 10<sup>-3</sup></td><td>1.2115</td></tr><tr><td>256</td><td>0.0027</td><td>1.6186</td><td>0.0115</td><td>1.6236</td><td>0.0198 × 10<sup>-3</sup></td><td>1.5780</td></tr><tr><td>512</td><td>0.0007</td><td>1.8754</td><td>0.0031</td><td>1.8772</td><td>0.0056 × 10<sup>-3</sup></td><td>1.8329</td></tr><tr><td>1024</td><td>0.0002</td><td>1.9785</td><td>0.0008</td><td>1.9820</td><td>0.0014 × 10<sup>-3</sup></td><td>1.9427</td></tr><tr><td rowspan="5">TrZ-scheme</td><td>64</td><td rowspan="5">h</td><td>0.0268</td><td>—</td><td rowspan="5">hu</td><td>0.1175</td><td>—</td><td rowspan="5">z</td><td>0.1824 × 10<sup>-3</sup></td><td>—</td></tr><tr><td>128</td><td>0.0083</td><td>1.6955</td><td>0.0354</td><td>1.7320</td><td>0.0550 × 10<sup>-3</sup></td><td>1.7296</td></tr><tr><td>256</td><td>0.0027</td><td>1.6182</td><td>0.0115</td><td>1.6248</td><td>0.0183 × 10<sup>-3</sup></td><td>1.5890</td></tr><tr><td>512</td><td>0.0007</td><td>1.8755</td><td>0.0031</td><td>1.8781</td><td>0.0050 × 10<sup>-3</sup></td><td>1.8734</td></tr><tr><td>1024</td><td>0.0002</td><td>1.9782</td><td>0.0008</td><td>1.9818</td><td>0.0013 × 10<sup>-3</sup></td><td>1.9792</td></tr></tbody></table>

namely the AcTrZ-scheme. On the other hand, the TrZ-method appears to be the less diffusive scheme but it is the only one which presents some oscillations in the solution. Of course, they could be related to the fact that the TrZ-method entails a decoupled numerical approximation of the shallow-water-Exner system, even if this kind of oscillations are not observed in the first-order version of TrZ-scheme. However, we underline that these oscillations reduce themselves when refining the mesh, refer to the previous work [13]. Finally, the AcZ-scheme is slightly more diffusive than the TrZ one but less than the AcTrZ-method.

## 7.3 Transient Riemann problems

Next, we consider three different transient Riemann problems to test the ability of our schemes to reproduce the correct solutions in this kind of situation as well. In table 2 we insert the initial conditions and the value for the coefficient $A_g$ present in the Grass formulation. We underline that in the third Riemann problem, $A_g$ is not constant anymore as it depends on the water height $h$. Moreover, we take $\rho_0 = 0.4$, $M = 200$ cells and $t_{end} = 0.2$s as ending time. We refer to the works [39, 32] for more details and for the analytical solutions reported in figures 3, 4, 5. In particular, we plot the free surface and bed elevation outputs using both first and second order method. Let us see the details. For the first RP (test A) represented in figure 3, on the left side we insert the AcTrZ solution while the TrZ one is shown on the right. We can generally observe that the numerical outputs follow closely the reference one even if the TrZ solution presents some small oscillations in correspondence of the middle discontinuity. The AcZ output is not inserted as it is very close to the AcTrZ one, even if less diffusive. Similar observations can be inferred for test B, thus we only show the AcTrZ solution on the left side of figure 4. We underline the neither the AcTrZ nor the AcZ solutions produce oscillations for these two RPs at both first and second order of accuracy. Finally, let us see test C outputs. On the right side of figure 4, we inserted the AcTrZ solution both at first and second order of accuracy. Once again, the numerical solution appear to
<figure><img src="image_8.png" /><figcaption>Figure 1: RP: dam break on movable bottom; free surface (left-top) and bed level (left-bottom and right). Comparison between the first (magenta) and second-order (blue) LP-AcZ scheme with $M = 100$ cells (left) and $M = 2000$ cells (right). Reference solution in red line.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 2: RP: dam break on movable bottom, bed level. Comparison among the three second-order schemes AcTrZ (magenta symbol), AcZ (yellow symbol) and TrZ (blue line) with $M = 100$ cells. Reference solution in red line.</figcaption></figure>
<table><caption>Table 2: Data for the transient Riemann problems.</caption><thead><tr><th>Test</th><th>h<sub>L</sub></th><th>h<sub>R</sub></th><th>u<sub>L</sub></th><th>u<sub>R</sub></th><th>z<sub>L</sub></th><th>z<sub>R</sub></th><th>A<sub>g</sub></th></tr></thead><tbody><tr><td>A</td><td>2</td><td>2</td><td>0.25495</td><td>2.3247449</td><td>3.0</td><td>2.846848</td><td>0.01</td></tr><tr><td>B</td><td>2.25</td><td>1.18868612</td><td>0.2050</td><td>2.4321238</td><td>5.0</td><td>5.124685</td><td>0.01</td></tr><tr><td>C</td><td>6</td><td>5.2</td><td>0.30037</td><td>15.16725</td><td>3.0</td><td>4.631165</td><td>0.01/h</td></tr></tbody></table>

reproduce correctly the reference one. However, while neither the first-order AcTrZ nor the first-order AcZ show any oscillations, the latter can be observed when using the second-order schemes. Thus, in figure 5 we compare the three second-order solutions and we zoom in the areas of interest to show that some small perturbations are present in the numerical outputs, probably due to the fact that less numerical diffusion is present. However, the correct solution is generally reproduced, even if $A_g$ is not constant anymore. We conclude saying that we verified that the numerical outputs converge to the analytical one when refining the mesh.

**Remark 2.** Moreover, even if without reporting here the data, even bigger values of $A_g$ have been considered to simulate highly erosive flow for the first dam break problem (section 7.2) and these three transient test cases. The interest in trying bigger values of $A_g$ resides in the fact that greater values means a change in the flow structure of the coupled model and consequently more instabilities could appear if the numerical scheme is decoupled. Using the value $A_g = 0.5$, we indeed verified the presence of great oscillations when using the TrZ method. On the other hand, the other two first-order methods AcTrZ and AcZ do not produce any oscillations and remain stable. However, some oscillations have been observed using their second-order extension, even if sufficiently controlled to keep the schemes stable. One could envisage a strategy to remove them by combining the first-order and second-order version in order to keep the second-order of accuracy in the stable parts of the solution, resembling a kind of flux-limiter approach (if a priori) or a MOOD approach (if a posteriori).

## 7.4 Sub-critical test case

For the following two numerical tests we refer to paper [20]. As initial condition we consider the sub-critical steady state

$$ \begin{cases} hu(x, t = 0) = 0.5 \\ z(x, t = 0) = 0.1(1 + e^{-(x-5)^2}) \\ \frac{u^2}{2} + g(h + z) = 6.386, \end{cases} $$

while the length of the channel is $L = 10.0m$. In figure 6 we insert the results for $A_g = 0.05$ and $A_g = 0.007$. In the latter case, AcZ-solution is not inserted for the sake of clarity as it is very similar to the AcTrZ-one. We observe that the three schemes give similar solutions, which confirm the observations of the previous test case; in order, from the least diffusive to the most diffusive, we have TrZ, AcZ and AcTrZ-scheme. It is also important to underline that, in these two numerical simulations,
<figure><img src="image_10.png" /><figcaption>Figure 3: Test A; free surface (up) and bed elevation (bottom) computed with AcTrZ method (left) and TrZ scheme (right). Reference solution (red line), first-order solution (blue symbol -.) and second-order solution (black dashed line). M = 200 cells.</figcaption></figure>

no oscillations appear even if, in the work [20], it has been explained that decoupled methods usually present oscillations in these two test cases.

## 7.5 "Lake at rest" solution and perturbation

Referring to [16], here we test the ability of the schemes to preserve the "lake at rest" steady state. First of all we consider the following stationary solution where $u = 0$, $h(x, t = 0) + z(x, t = 0) = 3m$ and

$$z(x, t = 0) = \begin{cases} 2 + 0.25(\cos(10\pi(x - 0.5)) + 1) & \text{if } 1.4 < x < 1.6 \\ 2 & \text{otherwise.} \end{cases}$$

The length of the channel is $L = 2.0m$. The numerical schemes were able to preserve this steady state up to an error of order $10^{-15}$. Then, we introduce some small perturbations in the initial data, namely we impose

$$h(x, t = 0) = \begin{cases} 3 - z(x, t = 0) + 0.001 & \text{if } 1.1 < x < 1.2 \\ 3 - z(x, t = 0) & \text{otherwise.} \end{cases}$$

In figure 7 we compare the results of first and second-order AcZ-scheme against the reference solution. We show the outputs of only the AcZ-method as the other two schemes give analogous results. We observe that the outcomes are indeed satisfying as they are in agreement with the ones showed in [16] and no unphysical oscillations appear.
<figure><img src="image_11.png" /><figcaption>Figure 4: Test B (left) and C (right); free surface (up) and bed elevation (bottom) computed with AcTrZ method. Reference solution (red line), first-order solution (blue pointed line) and second-order solution (black dashed line). M = 200 cells.</figcaption></figure>

## 7.6 "Constant bed slope" equilibrium for steady flow regimes

In this section we numerically show that our method is capable of preserving the "constant bed slope" equilibrium and that the solution evolves to the steady state if steady boundary conditions are imposed. We refer to [32] and we take $\mu_f = 0.020006460818026 \text{ s m}^{-1/3}$, $A_g = 0.01$, $L = 100 \text{ m}$ and $M = 100$ cells so that $\Delta x = 1 \text{ m}$. Then, for the equilibrium, we consider the following slope $S_{eq} = -0.002$ for $z$. Our numerical methods are able to preserve it with an error machine of $10^{-12}$. We underline that at the left and right boundaries, for the variable $z$ we need to impose that the slope is constant. For instance, for the right boundary, this means that the value in the ghost cell is given by $z_{\text{out}} - z_{\text{end}} = z_{\text{end}} - z_{\text{end}-1}$, where by $z_{\text{end}}$ and $z_{\text{end}-1}$ we indicate the value of $z$ in the last and second-last cell.

Then, we move away from the steady state and consider as initial condition $h(x, t=0) = 0.943 \text{ m}$, $q(x, t=0) = 1 \text{ m}^3/\text{s}$ and either $S_0 = -0.007$ or $S_0 = 0$. As for the boundary conditions, at the inlet we impose $z_{\text{In}} = 2 \text{ m}$ and $q_{\text{In}} = 1 \text{ m}^3/\text{s}$, while at the outlet we use $h_{\text{Out}} = 0.943 \text{ m}$. Once again, at the right boundary, for the variable $z$ we imposed that the slope is constant. Then, in figure 8 we insert the results using $S_0 = -0.007$ (left) and $S_0 = 0$ (right). In both cases, we observe that the numerical solution converges towards the exact one, namely the "constant bed slope" equilibrium. The outputs have been computed with the first-order LP-AcTrZ scheme, but all the numerical methods we presented give analogous results.
<figure><img src="image_12.png" /><figcaption>Figure 5: Test C; free surface (up) and bed elevation (bottom) computed with the three second order schemes. AcTrZ (blue line), AcZ (black dashed line), TrZ (magenta pointed line) and reference solution (red line). M = 200 cells.</figcaption></figure>

## 7.7 Dam break with experimental values

Here we present the last numerical test, in which we compare our numerical solution against experimental data. We consider once again a Riemann problem performed at the Université catholique de Louvain with initial condition given by zero-velocity, $h_L = 0.1$ m, $h_R = 1e-3$ m for the water height and flat topography. Then we take $L = 2.5$ m, MPM formula (1.3) for the solid transport discharge with $d = 3.2$ mm, $\rho_0 = 0.4$, $s = 0.540$, $\mu_f = 0.03$ and $\theta_c^* = 0.045$. Refer to [27] for more details about this experiment. In figure 9 we insert the free surface and bed elevation numerical and experimental outputs computed at different times $t_{end} = 5t_0, 7.5t_0, 10t_0$ s with $t_0 = \sqrt{gh_0} \approx 0.101$. In particular, we used the AcZ scheme but analogous solution can be found using the AcTrZ or TrZ method, where once again the TrZ solution would be the less diffusive while the AcTrZ output the most diffusive. The results are considered satisfying as they appear to match the experimental data and are comparable to the ones obtained in [27]. Note that the difference between the numerical and experimental output could be related to the fact that we are both neglecting erosion processes (meaning that smaller fractions of sediment could be in suspension into the water) and non-hydrostatic effects, refer also to [25]. Then, when using the second-order schemes we observed the production of some small oscillations which were dumped as time went on. Thus, in figure 11 we inserted the second-order solutions using the AcTrZ and AcZ schemes at times $t_{end} = 5t_0, 7.5t_0, 10t_0$ s, while in figure 10 at time $t_{end} = 1.5t_0$ to better show the presence of oscillations. Note that we did not inserted the second-order TrZ solution as it gave complex number. Concluding, even if small oscillations are present, probably due to the fact that at second-order of accuracy there is less diffusion, the solutions are considered satisfying. Moreover, at this stage we
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Figure 6: Flow over a movable bump; free surface and bed level. Values $A_g = 0.005$ (up) and $A_g = 0.07$ (bottom). Schemes AcTrZ (magenta), AcZ (yellow) and TrZ (blue) with $M = 100$ cells.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 7: Propagation of perturbation; bed level $z$ and free surface $z + h$ (top), velocity (bottom). Comparison between the first (magenta) and second-order (blue) LP-AcZ scheme with $M = 200$ cells. Reference solution in red line.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 8: "Constant bed slope" equilibrium for steady flow regimes; bed level $z$ and free surface $z + h$. Initial condition (blue line), numerical steady state (green symbol) and exact steady state (red line). LP-AcTrZ scheme with $M = 100$ cells.</figcaption></figure>
are still neglecting source terms related to the erosion phenomena of the sediment, which could help stabilizing the numerical output, see [25].

# 8 Conclusions

In the previous work [13], it has been introduced a well-balanced second-order Lagrange-projection scheme for the shallow water Exner model which takes into account the topography $z$ only in the projection step. In the present work, we described two more different ways of discretizing the bed level, by considering $z$ only in the acoustic step or in both of them. Both methods have been constructed in such a way to be well-balanced and second-order accurate at the same time. Here by well-balanced we mean that the schemes preserve the "lake at rest" and the "constant bed slope" equilibrium steady states. Moreover, both Grass and MPM formulation have been considered. We tested these methods considering several numerical experiments and we generally observed that when using their first-order version no oscillations are present in the numerical results. When moving to second-order of accuracy however this is not always true. In particular, the TrZ scheme ($z$ only in the transport step) revealed itself to be the worst method among the three in this sense. On the other hand, the other two methods AcTrZ and AcZ performed a lot better, showing only minor and controlled oscillations in some of the numerical tests. Further improvements could be related to the extension of these numerical schemes to two dimensions (following the lines of [13]), higher order of accuracy and the use of other formulations for the solid transport discharge. Finally, it would be interesting to consider the implicit version of such numerical methods, especially in situations of weak interaction between the flow and the sediments, where long-time simulations are needed. Indeed, such a topic is already subject of an ongoing work.

# Acknowledgments

This work was partially supported by a grant from Région Île-de-France.

We also would like to thank Professor M. J. Castro and Professor T. Morales for the valuable discussions and for the experimental data used in section 7.7.

# References

[1] Emmanuel Audusse, Christophe Berthon, Christophe Chalons, Olivier Delestre, Nicole Goutal, Magali Jodeau, Jacques Sainte-Marie, Jan Giesselmann, Georges Sadaka. *Sediment transport modelling: Relaxation schemes for Saint-Venant – Exner and three layer models*. ESAIM: Proc. 38 78-98 (2012). DOI: 10.1051/proc/201238005

[2] E. Audusse, F. Bouchut, M.O. Bristeau, R. Klein and B. Perthame. *A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows*. SIAM Journal on Scientific Computing, 25:2050–2065, 2004.
<figure><img src="image_17.png" /><figcaption>Figure 9: Dam break with experimental values; numerical bed level *z* (red line), numerical free surface *z* + *h* (blue line) and experimental data (red symbol). Solution at different times: *t* = 5*t*₀ *s* (up), *t* = 7.5*t*₀ *s* (middle), *t* = 10*t*₀ *s* (bottom). LP-AcZ scheme with *M* = 200 cells.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 10: Dam break with experimental values; AcTrZ scheme (blue line) and AcZ scheme (black dashed line). Solution at time $t = 1.5t_0$ s. $M = 200$ cells.</figcaption></figure>

[3] E. Audusse, C. Chalons, P. Ung. *A simple three-wave approximate Riemann solver for the Saint-Venant-Exner equations*. Int J Numer Meth Fluids. 2018; 87:508–528.

[4] Christophe Berthon, Stéphane Cordier, Minh Le, Olivier Delestre. *An analytical solution of Shallow Water system coupled to Exner equation*. Comptes Rendus Mathématique, Elsevier Masson, 2012, 350 (3-4), pp.183-186.

[5] François Bouchut, Christophe Chalons, Sébastien Guisset. *An entropy satisfying two-speed relaxation system for the barotropic Euler equations. Application to the numerical approximation of low Mach number flows*. Numerische Mathematik, Springer Verlag, 2020, 145, pp.35-76.
ff10.1007/s00211-020-01111-5ff. ffhal-01661275v2f

[6] R. Briganti, N. Dodd, D. Kelly, and D. Pokrajac. (2012). *An efficient and flexible solver for the simulation of the morphodynamics of fast evolving flows on coarse sediment beaches*. International Journal for Numerical Methods in Fluids. 69. 859 - 877. 10.1002/fld.2618.

[7] M. J. Castro Díaz, E. D. Fernández-Nieto and A.M. Ferreiro. *Sediment transport models in Shallow Water equations and numerical approach by high order finite volume methods*. Computers & Fluids, Volume 37, Issue 3, 2008, Pages 299-316.

[8] Christophe Berthon, Françoise Foucher. *Efficient well-balanced hydrostatic upwind schemes for shallow-water equations*. 2012. Journal of Computational Physics. 231. 4993–5015. 10.1016/j.jcp.2012.02.031.

[9] Marie Billaud Friess, Benjamin Boutin, Filipa Caetano, Gloria Faccaroni, Samuel Kokh, et al.. *A second order anti-diffusive Lagrange-remap scheme for two-component flows*. ESAIM: Proceedings, EDP Sciences, 2011, 32, pp.149-162. ff10.1051/proc/2011018ff. ffhal-00708605f
<figure><img src="image_19.png" /><figcaption>Figure 11: Dam break with experimental values; AcTrZ scheme (blue line), AcZ scheme (black dashed line) and experimental data (red symbol). Solution at different times: $t = 5t_0$ s (up), $t = 7.5t_0$ s (middle), $t = 10t_0$ s (bottom). M = 200 cells.</figcaption></figure>
[10] François Bouchut. *Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources*. Frontiers in mathematics, 2004.

[11] Manuel J. Castro Díaz, Christophe Chalons and Tomás Morales De Luna. *A fully well-balanced Lagrange-Projection type scheme for the Shallow-water equations*. 2018. SIAM J. Numer. Anal., 56(5), 3071–3098.

[12] M.J. Castro, T. Morales de Luna and C. Parés. *Well-Balanced Schemes and Path-Conservative Numerical Methods*. Handbook of Numerical Analysis, Vol. 18. 2017. Pages 131-175.

[13] C. Chalons and A. Del Grosso. *A second-order well-balanced numerical scheme for Shallow Water Exner equations in 1D and 2D*.

[14] Christophe Chalons, Mathieu Girardin, Samuel Kokh. *An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes*. Journal of Computational Physics, Elsevier, 2017, 335, pp.885-904. 10.1016/

[15] Christophe Chalons, Mathieu Girardin, Samuel Kokh. (2014). *An All-Regime Lagrange-Projection Like Scheme for the Gas Dynamics Equations on Unstructured Meshes*. Communications in Computational Physics. 20. 10.4208/cicp.260614.061115a.

[16] Christophe Chalons, Pierre Kestener, Samuel Kokh, and Maxime Stauffert. *A large time-step and well-balanced Lagrange-Projection type scheme for the Shallow-water equations*. 2016. Communications in Mathematical Sciences. 15. 10.4310/CMS.2017.v15.n3.a9.

[17] Christophe Chalons, Samuel Kokh, Mathieu Girardin. (2013). *Large Time Step and Asymptotic Pre-serving Numerical Schemes for the Gas Dynamics Equations with Source Terms*. SIAM Journal on Scientific Computing. 35. 10.1137/130908671.

[18] F. Coquel, E. Godlewski, B. Perthame, A. In, and P. Rascle. *Some new Godunov and relaxation methods for two-phase flow problems*. Godunov methods (Oxford, 1999), pages 179–188, 2001.

[19] F. Coquel and B. Perthame. *Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics*. SIAM Journal on Numerical Analysis, 35(6):2223–2249, 1998.

[20] Stéphane Cordier, Minh Le, Tomas Morales de Luna. *Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help*. Advances in Water Resources, 34(8): 980–989. August 2011.

[21] Frédéric Duboc, Cédric Enaux, Stéphane Jaouen, Hervé Jourdren, Marc Wolff. *High-order dimensionally split Lagrange-remap schemes for compressible hydrodynamics*. Comptes Rendus Mathematique Volume 348, Issues 1–2, January 2010, Pages 105-110.

[22] F. Franzini, S. Soares-Frazão. Coupled finite-volume scheme with adapted Augmented Roe scheme for simulating morphological evolution in arbitrary cross-sections. J. Hydroinform. 2018, 20, 1111–1130.
[23] G. Gallice. *Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source.*
C. R. Math. Acad. Sci. Paris 334, no. 8, pp. 713-716, (2002).

[24] G. Gallice. *Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations
in Lagrangian or Eulerian coordinates*. Numer. Math. 94 no. 4, pp. 673-713, (2003).

[25] J. C. González-Aguirre, M.J. Castro, T. Morales de Luna. *A robust model for rapidly varying flows
over movable bottom with suspended and bedload transport: modelling and numerical approach.*
*Advances in Water Resources* Volume 140, June 2020, 103575.

[26] Sigal Gottlieb and Chi-Wang Shu. *Total variation diminishing RUNGE-KUTTA schemes*. 1996. Math-
ematics of Computation. 67. 10.1090/S0025-5718-98-00913-2.

[27] Harry Putu Gunawan. *Numerical simulation of shallow water equations and related models*. Gen-
eral Mathematics [math.GM]. Université Paris-Est, 2015. English. ffNNT : 2015PEST1010ff. ff-
tel01216642v2f

[28] J. Hudson. *Numerical techniques for morphodynamic modelling*. PhD thesis, University of Reading,
October 2001.

[29] R.J. LeVeque. *Numerical Methods for Conservation Laws*. 1992 - Springer.

[30] Xin Liu, Abdolmajid Mohammadian, Alexander Kurganov, Julio Angel Infante Sedano. *Well-balanced central-upwind scheme for a fully coupled shallow water system modeling flows over erodible bed*. Journal of Computational Physics, Volume 300, 1 November 2015, Pages 202-218.

[31] S. Martínez-Aranda, J. Murillo, P. García-Navarro. *A comparative analysis of capacity and non-capacity formulations for the simulation of unsteady flows over finite-depth erodible beds*. Advances in Water Re-sources 130: 91-112 (2019).

[32] S. Martínez-Aranda, J. Murillo, P. García-Navarro. *A 1D numerical model for the simulation of unsteady and highly erosive flows in rivers*. Computers & Fluids 181, 8-34 (2019)

[33] S. Martínez-Aranda, J. Murillo, P. García-Navarro. *Comparison of new efficient 2D models for the simulation of bedload transport using the augmented Roe approach*. Advances in Water Resources 153, 103931 (2021).

[34] R. Meurice, S. Soares-Frazão. *A 2D HLL-based weakly coupled model for transient flows on mobile beds*. J. Hydroinform. 2020, 22, 1351-1369.

[35] Victor Michel-Dansac, Christophe Berthon, Stéphane Clain, Françoise Foucher. *A well-balanced scheme for the shallow-water equations with topography*. Computers and Mathematics with Appli-cations. 72. 10.1016/j.camwa.2016.05.015. (2015)

[36] Victor Michel-Dansac, Christophe Berthon, Stéphane Clain, Françoise Foucher. *A well-balanced scheme for the shallow-water equations with topography or Manning friction*. Journal of Computa-tional Physics, Elsevier, 2017, 335, pp.115-154.
[37] Tomás Morales De Luna, Manuel J. Castro Díaz and Christophe Chalons. *High order fully well-balanced Lagrange-Projection scheme for Shallow-water*. Commun. Math. Sci., Vol. 18, No. 3, pp. 781–807 (2020)

[38] Tomás Morales de Luna, Manuel J. Castro Díaz, Carlos Parés Madroñal, Enrique D. Fernández-Nieto. *On a Shallow Water Model for the Simulation of Turbidity Currents*. (2009). Communications in Computational Physics - COMMUN COMPUT PHYS. 6. 10.4208/cicp.2009.v6.p848.

[39] J. Murillo and P. García-Navarro. *An Exner-based coupled model for two-dimensional transient flow over erodible bed*. Journal of Computational Physics, Volume 229, Issue 23, p. 8704-8732. Nov 2010.

[40] Enrique Fernández-Nieto, Carine Lucas, Tomás Morales de Luna, Stephane Cordier. On the influence of the thickness of the sediment moving layer in the definition of the bedload transport formula in Exner systems. Computers & Fluids 91, 87–106 (2014).

[41] E. F. Toro. *Riemann Solvers and Numerical Methods for Fluid Dynamics*, Third Edition. Springer-Verlag, 2009.

[42] C. B. Vreugdenhil. *Numerical Methods for Shallow Water Flow*. Dordrecht ; Boston: Kluwer Academic Publishers (1994). Part of the Water Science and Technology Library book series (WSTL, volume 13).