# A fully well-balanced lagrange-projection-type scheme for the shallow-water equations

Manuel Jesús Castro Díiaz, Christophe Chalons, Tomás Morales de Luna

► To cite this version:

Manuel Jesús Castro Díiaz, Christophe Chalons, Tomás Morales de Luna. A fully well-balanced lagrange-projection-type scheme for the shallow-water equations. SIAM Journal on Numerical Analysis, 2018, 56 (5), pp.3071-3098. 10.1137/17M1156101 . hal-03993092

HAL Id: hal-03993092

https://hal.science/hal-03993092

Submitted on 28 Feb 2023

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
A FULLY WELL-BALANCED LAGRANGE-PROJECTION TYPE
SCHEME FOR THE SHALLOW-WATER EQUATIONS*

MANUEL J. CASTRO DÍAZ†, CHRISTOPHE CHALONS‡, AND TOMÁS MORALES DE LUNA§

**Abstract.** This work focuses on the numerical approximation of the Shallow Water Equations (SWE) using a Lagrange-Projection type approach. We propose a fully well-balanced explicit and positive scheme using relevant reconstruction operators. By fully well-balanced, it is meant that the scheme is able to preserve stationary smooth solutions of the model with non zero velocity, including of course also the well-known lake at rest equilibrium. Numerous numerical experiments illustrate the good behaviour of the scheme.

**Key words.** Shallow water equations, Finite Volume Method, Lagrange-Projection scheme

**AMS subject classifications.** 74S10, 35L60, 35L65, 74G15

**1. Introduction.** We are interested in the numerical approximation of the well-known Shallow Water Equations (SWE), given by

$$ (1) \qquad \left\{ \begin{aligned} \partial_t h + \partial_x (hu) &= 0, \\ \partial_t (hu) + \partial_x \left(hu^2 + g \frac{h^2}{2}\right) &= -gh \partial_x z, \end{aligned} \right. $$

where $z(x)$ denotes a given smooth topography and $g > 0$ is the gravity constant. The primitive variables are the water depth $h \ge 0$ and its velocity $u$, which both depend on the space and time variables, respectively $x \in \mathbb{R}$ and $t \in [0, \infty)$. At time $t = 0$, we assume that the initial water depth $h(x, t = 0) = h_0(x)$ and velocity $u(x, t = 0) = u_0(x)$ are given. In order to shorten the notations, we will use the following condensed form of (1), namely

$$ (2) \qquad \partial_t \mathbf{U} + \partial_x \mathbf{F}(\mathbf{U}) = \mathbf{S}(\mathbf{U}, z), $$

where $\mathbf{U} = (h, hu)^T$, $\mathbf{F}(\mathbf{U}) = (hu, hu^2 + gh^2/2)^T$ and $\mathbf{S}(\mathbf{U}, z) = (0, -gh\partial_x z)^T$. In case of necessity and since the topography is assumed to be constant in time in this work, we will also make use of the following equivalent form of (1), namely

$$ (3) \qquad \left\{ \begin{aligned} \partial_t h + \partial_x (hu) &= 0, \\ \partial_t (hu) + \partial_x \left(hu^2 + g \frac{h^2}{2}\right) + gh \partial_x z &= 0, \\ \partial_t z &= 0, \end{aligned} \right. $$

<footnote>*Submitted to the editors Nov 2017.</footnote>

<footnote>**Funding:** This research has been partially supported by the Spanish Government and FEDER through the coordinated Research projects MTM 2015-70490-C2-1-R and MTM 2015-70490-C2-2-R. The authors were partially supported by the Innovative Training Networks (ITN) grant 642768 (Mod-CompShock). The second author was partially supported by the ANR project ACHYLLES (ANR-14-CE25-0001-03).</footnote>

<footnote>†Departamento de Análisis Matemático, Facultad de Ciencias, Universidad de Málaga, 29071 Málaga, Spain (castro@anamat.cie.uma.es).</footnote>

<footnote>‡Laboratoire de Mathématiques de Versailles, UVSQ, CNRS, Université Paris-Saclay, 78035 Versailles, France (christophe.chalons@uvsq.fr).</footnote>

<footnote>§Dpto. de Matemáticas. Universidad de Córdoba. Campus de Rabanales. 14071 Córdoba, Spain (tomas.morales@uco.es).</footnote>
and $\mathbf{V} = (h, hu, z)^T$. At last, recall that (1) is strictly hyperbolic over the phase space $\Omega = \{(h, hu)^T \in \mathbb{R}^2 | h > 0\}$ and that the eigenstructure is composed by two genuinely non linear characteristic fields associated with the eigenvalues $\{u - c, u + c\}$ where the sound speed is defined by $c = \sqrt{gh}$. We recall also that the regions where $u^2 < c^2$ (resp. $u^2 > c^2$) are called subcritical (resp. supercritical).

As already stated, we are especially interested in the design of a numerical scheme satisfying the so-called fully well-balanced property, meaning that it should be able to preserve discrete approximations of the smooth stationary solutions of (1). These steady states are governed by the ordinary differential system $\partial_x F(\mathbf{U}) = \mathbf{S}(\mathbf{U}, z)$, namely

$$ (4) \qquad hu = \text{constant}, \quad \frac{u^2}{2} + g(h+z) = \text{constant}. $$

Recall that the well-known "lake at rest" solution corresponds to

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

and that numerical schemes preserving these particular stationary solutions are said to be well-balanced (see for instance the recent book [26] for a review). Here we are also interested in preserving equilibriums (4) with a non-zero velocity, leading to the so-called fully well-balanced property.

In order to reach this objective, we propose to use a strategy based on a Lagrange-Projection (or equivalently Lagrange-Remap) decomposition, which allows to naturally decouple the acoustic and transport terms of the model. Such a decomposition proved to be useful and very efficient to deal with subsonic or near low-Froude number flows. In this case, the usual CFL time step limitation of Godunov-type schemes is indeed driven by the acoustic waves and can thus be very restrictive. The Lagrange-Projection strategy allows to design a very natural implicit-explicit and large time-step schemes with a CFL restriction based on the (slow) transport waves and not on the (fast) acoustic waves. In this paper, we will stay focused on explicit schemes but the reader is referred to the pioneering work [20] and the more recent ones [14], [15], [16], [17] for more details. Note that the large time step Lagrange-Projection scheme proposed in [17] for the SWE is well-balanced, but not fully well-balanced.

There is a huge amount of works about the design of numerical schemes for the SWE, but most of them intend to satisfy the well-balanced property, but not the full well-balanced property. To mention only a few of them, we refer for instance the reader to the following well-known contributions [3], [27], [30], [22], [23], [1], [19], [6], [7], [2]. We also refer to the books [9] and [26] as well as [29] and [13] which provide additional references and overviews. As far as the full well-balanced property is concerned, let us quote for instance, [25], [34], [36], [11], [10], [35], [33], [4], [5], [12], [28], [31].

The outline of the paper is as follows. In Section 2 we give a general introduction of the Lagrange-Projection approach in the case of flat topography. This technique will be generalized in Section 3 for the case of arbitrary topography in order to develop a fully well-balanced scheme. Finally, in Section 4, some numerical results will be shown in order to test this new scheme.

## 2. The Lagrange-Projection decomposition.
The aim of this section is to briefly present the Lagrange-Projection splitting strategy considered in this paper and to introduce the basics of the numerical approximation.

Let us begin with the Lagrange-Projection decomposition. By using the chain rule for the space derivatives, one first writes system (1) under the following equivalent form
72 for smooth solutions, namely

$$73 \quad \left\{ \begin{array}{l} \partial_t h + u \partial_x h + h \partial_x u = 0 \\ \partial_t (hu) + u \partial_x (hu) + hu \partial_x u + \partial_x (g \frac{h^2}{2}) = -gh \partial_x z. \end{array} \right. $$

74 Then, the splitting simply consists in first solving the Lagrangian form of this system, 75 namely

$$76 \quad (6) \quad \left\{ \begin{array}{l} \partial_t h + h \partial_x u = 0, \\ \partial_t (hu) + hu \partial_x u + \partial_x (g \frac{h^2}{2}) = -gh \partial_x z, \end{array} \right. $$

which gives in Lagrangian coordinates, with $\tau = 1/h$ and $\tau \partial_x = \partial_m$,

$$77 \quad \left\{ \begin{array}{l} \partial_t \tau - \partial_m u = 0, \\ \partial_t u + \partial_m (g \frac{h^2}{2}) = -gh \partial_m z, \end{array} \right. $$

78 and then the transport system

$$79 \quad (7) \quad \left\{ \begin{array}{l} \partial_t h + u \partial_x h = 0, \\ \partial_t (hu) + u \partial_x (hu) = 0. \end{array} \right. $$

80 *Remark 2.1.* This splitting procedure can be interpreted as a two-step algorithm 81 consisting in solving first the system in Lagrangian coordinates and then projecting the 82 results in the Eulerian coordinates. See [24] for further details.

83 **2.1. The Lagrange-Projection numerical algorithm.** Space and time will be 84 discretized using a space step $\Delta x$ and a time step $\Delta t$ into a set of cells $[x_{j-1/2}, x_{j+1/2})$ 85 and instants $t^{n+1} = n\Delta t$, where $x_{j+1/2} = j\Delta x$ and $x_j = (x_{j-1/2} + x_{j+1/2})/2$ are 86 respectively the cell interfaces and cell centers, for $j \in \mathbb{Z}$ and $n \in \mathbb{N}$. For a given initial 87 condition $x \mapsto \mathbf{U}^0(x)$, we will consider a discrete initial data $\mathbf{U}_j^0$ which approximates 88 $\frac{1}{\Delta x} \int_{x_{j-1/2}}^{x_{j+1/2}} \mathbf{U}^0(x) dx$, for $j \in \mathbb{Z}$. Therefore, the proposed algorithm aims at computing 89 an approximation $\mathbf{U}_j^n$ of $\frac{1}{\Delta x} \int_{x_{j-1/2}}^{x_{j+1/2}} \mathbf{U}(x, t^n) dx$ where $x \to \mathbf{U}(x, t^n)$ is the exact solution 90 of (1) at all time $t^n$, $n \in \mathbb{N}$. Given the sequence $\{\mathbf{U}_j^n\}_j$, it is a matter of defining the 91 sequence $\{\mathbf{U}_j^{n+1}\}_j$, $n \in \mathbb{N}$, since $\{\mathbf{U}_j^0\}_j$ is assumed to be known.

92 Using these notations, the overall Lagrange-Projection algorithm can be described 93 as follows: for a given discrete state $\mathbf{U}_j^n = (h, hu)_j^n$, $j \in \mathbb{Z}$ that describes the system at 94 instant $t^n$, the computation of the approximation $\mathbf{U}_j^{n+1} = (h, hu)_j^{n+1}$ at the next time 95 level is a two-step process defined by

96 1. Update $\mathbf{U}_j^n$ to $\mathbf{U}_j^{n+1-}$ by approximating the solution of (6),

97 2. Update $\mathbf{U}_j^{n+1-}$ to $\mathbf{U}_j^{n+1}$ by approximating the solution of (7).

98 **2.2. The case of a flat topography.** In order to prepare the forthcoming devel- 99 opments, we first consider the simple situation of a flat topography, which corresponds 100 to the case $\partial_x z = 0$. In this case and using standard notations, the Lagrangian step 101 writes (see [24])

$$102 \quad \left\{ \begin{array}{l} \tau_j^{n+1-} = \tau_j^n + \frac{\Delta t}{h_j^n \Delta x} (u_{j+1/2}^* - u_{j-1/2}^*), \\ u_j^{n+1-} = u_j^n - \frac{\Delta t}{h_j^n \Delta x} (p_{j+1/2}^* - p_{j-1/2}^*), \end{array} \right. $$
with

$$u_{j+1/2}^* = u^*(\mathbf{U}_j^n, \mathbf{U}_{j+1}^n) := \frac{1}{2}(u_j^n + u_{j+1}^n) - \frac{1}{2a_{j+1/2}}(p_{j+1}^n - p_j^n),$$

$$p_{j+1/2}^* = p^*(\mathbf{U}_j^n, \mathbf{U}_{j+1}^n) := \frac{1}{2}(p_j^n + p_{j+1}^n) - \frac{a_{j+1/2}}{2}(u_{j+1}^n - u_j^n),$$

and $p_j^n = g(h_j^n)^2/2$ for all $j$, or equivalently

$$ (8) \qquad \left\{ \begin{array}{l} L_j^n h_j^{n+1-} = h_j^n, \\ L_j^n (hu)_j^{n+1-} = (hu)_j^n - \frac{\Delta t}{\Delta x} (p_{j+1/2}^* - p_{j-1/2}^*), \end{array} \right. $$

with

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

Here, the constant $a_{j+1/2}$ has to be chosen sufficiently large for the sake of stability, and more precisely larger than the Lagrangian sound speed $hc$ according to the well-known subcharacteristic condition. In practice it is required that $a_{i+1/2}$ is greater than the values $h\sqrt{gh}$ at the interface. The definition of $a_{j+1/2}$ will be given later on and we also refer for instance the reader to [9, 24, 14, 16] for more details. Moreover, the Lagrangian step is stable under the Courant-Friedrichs-Lewy (CFL) condition

$$\Delta t \max_j \{\max(\tau_j^n, \tau_{j+1}^n) a_{j+1/2}\} \leq \frac{1}{2} \Delta x.$$

As far as the transport step is concerned, an upwind scheme is applied where the value $u_{j+1/2}^*$ is the velocity at the interface $x_{j+1/2}$.

$$ \left\{ \begin{array}{l} h_j^{n+1} = h_j^{n+1-} - \frac{\Delta t}{\Delta x} (u_{j+1/2}^{*-} (h_j^{n+1-} - h_j^{n+1-}) + u_{j-1/2}^{*+} (h_j^{n+1-} - h_{j-1}^{n+1-})), \\ (hu)_j^{n+1} = (hu)_j^{n+1-} - \frac{\Delta t}{\Delta x} (u_{j+1/2}^{*-} ((hu)_j^{n+1-} - (hu)_j^{n+1-}) + u_{j-1/2}^{*+} ((hu)_j^{n+1-} - (hu)_j^{n+1-})), \end{array} \right. $$

$$ \begin{array}{l} \text{where } u_{j\pm 1/2}^{*,+} = \max(u_{j\pm 1/2}^*, 0) \text{ and } u_{j\pm 1/2}^{*,-} = \min(u_{j\pm 1/2}^*, 0) \text{ for all } j. \\ \text{The CFL condition associated with the transport step reads} \end{array} $$

$$\Delta t \max_j \{(u_{j-1/2})^+ - (u_{j+1/2})^-\} \leq \Delta x.$$

Easy calculations show that the system can be written in the equivalent form

$$ \left\{ \begin{array}{l} h_j^{n+1} = L_j^n h_j^{n+1-} - \frac{\Delta t}{\Delta x} (h_{j+1/2}^{n+1-} u_{j+1/2}^* - h_{j-1/2}^{n+1-} u_{j-1/2}^*), \\ (hu)_j^{n+1} = L_j^n (hu)_j^{n+1-} - \frac{\Delta t}{\Delta x} ((hu)_{j+1/2}^{n+1-} u_{j+1/2}^* - (hu)_{j-1/2}^{n+1-} u_{j-1/2}^*), \end{array} \right. $$

where

$$h_{j+1/2}^{n+1-} = \begin{cases} h_j^{n+1-} & \text{if } u_{j+1/2}^* \ge 0, \\ h_{j+1}^{n+1-} & \text{if } u_{j+1/2}^* \le 0, \end{cases} \quad (hu)_{j+1/2}^{n+1-} = \begin{cases} (hu)_j^{n+1-} & \text{if } u_{j+1/2}^* \ge 0, \\ (hu)_{j+1}^{n+1-} & \text{if } u_{j+1/2}^* \le 0, \end{cases}$$

which gives that the whole scheme is conservative in the usual sense of finite volume methods and writes

$$ \left\{ \begin{array}{l} h_j^{n+1} = h_j^n - \frac{\Delta t}{\Delta x} (h_{j+1/2}^{n+1-} u_{j+1/2}^* - h_{j-1/2}^{n+1-} u_{j-1/2}^*), \\ (hu)_j^{n+1} = (hu)_j^n - \frac{\Delta t}{\Delta x} ((hu)_{j+1/2}^{n+1-} u_{j+1/2}^* + p_{j+1/2}^* - (hu)_{j-1/2}^{n+1-} u_{j-1/2}^* - p_{j-1/2}^*). \end{array} \right. $$
These formulas will be useful in the next sections to deal with a generalization of this scheme to the case of a non flat topography and impose the fully well-balanced property.

**2.3. An equivalent formulation.** In this section, we propose an equivalent formulation of the previous scheme based on the Reynolds' transport theorem which can be expressed for any scalar value $X$ as

$$ \frac{d}{dt} \int_{\Omega(t)} X(x,t) dV = \int_{\Omega(t)} \left( \frac{\partial}{\partial t} X(x,t) + \nabla \cdot (X(x,t)v(x,t)) \right) dV, $$

where $dV$ is a volume element at point $x$ (the integration variable), and $v(x,t)$ is the velocity at which the domain $\Omega(t)$ is moving. In the one dimension case considered in this paper, and considering that $\Omega(t)$ is an interval with limits $a(t)$ and $b(t)$, the Reynolds's theorem reduces to

$$ \frac{d}{dt} \int_{a(t)}^{b(t)} X(x,t) dx = \int_{a(t)}^{b(t)} \left( \frac{\partial}{\partial t} X(x,t) + \partial_x (X(x,t)v(x,t)) \right) dx. $$

This theorem expresses the time variation of the amount of $X$ contained inside the volume $\Omega(t)$ which moves at velocity $v(x,t)$.

*The Lagrangian coordinates.* Considering that the velocity $v(x,t)$ equals the fluid velocity $u(x,t)$ in the SWE, the Reynolds' theorem allows to describe the fluid motion in the Lagrangian coordinates, assuming that the observer moves with the fluid flow. Considering now successively that $X = 1, h, hu$ and denoting by $V(t)$ the infinitesimal control volume, we get

$$ (9) \qquad \frac{d}{dt} V(t) = V(t) \partial_x u(x(t; \xi), t), $$

$$ (10) \qquad \frac{d}{dt} \left( V(t)h(x(t; \xi), t) \right) = 0, $$ $$ \frac{d}{dt} \left( V(t)hu(x(t; \xi), t) \right) = -V(t) \partial_x p(x(t; \xi), t), $$

where each point $x(t; \xi)$ of the control volume is given by the solution of the following ordinary differential equation, with initial condition $\xi$,

$$ \left\{ \begin{array}{l} \partial_t x(t; \xi) = u(x(t; \xi), t) \\ x(0; \xi) = \xi. \end{array} \right. $$

Note that with the above notations, and $\xi_1 = a(0)$ and $\xi_2 = b(0)$,

$$ V(t) = \int_{x(t; \xi_1)}^{x(t; \xi_2)} 1 dx = \int_{\xi_1}^{\xi_2} \partial_{\xi} x(t; \xi) d\xi $$

so that for an infinitesimal control volume one can write

$$ V(t) = V(0) \partial_{\xi} x(t; \xi). $$

If we now denote $\bar{\alpha}(\xi, t) = \alpha(x(t; \xi), t)$ for any scalar quantity $\alpha$, the previous relations also write

$$ \frac{d}{dt} V(t) = V(0) \partial_{\xi} \bar{u}(\xi, t), $$
$$ \frac{d}{dt} (V(t) \bar{h}(\xi, t)) = 0, $$

114

$$ (11) \qquad \frac{d}{dt} \left( V(t) \overline{hu}(\xi, t) \right) = -V(0) \partial_{\xi} \bar{p}(\xi, t). $$

115

The Lagrangian step. Let us now turn back to the numerical setting and consider for a given $j$ that the control volume is defined at time $t = 0$ by the cell $(x_{j-1/2}, x_{j+1/2})$. The endpoints $x_{\pm}(t; x_{j\pm 1/2})$ of the control volume at a given time $t$ are thus defined by the solutions of the following ordinary differential equations

117

118

119

$$ (12) \qquad \left\{ \begin{array}{l} \partial_t x_{\pm}(t; x_{j\pm 1/2}) = u(x_{\pm}(t; x_{j\pm 1/2}), t) \\ x_{\pm}(0; x_{j\pm 1/2}) = x_{j\pm 1/2}. \end{array} \right. $$

120

Denoting $h_j^{n+1-}$ and $(hu)_j^{n+1-}$ constant approximate values of $h$ and $(hu)$ at time $t^{n+1}$ on the moving domain whose volume is $V_j(\Delta t)$, a simple way to discretize (9) would be

121

122

$$ (13) \qquad \left\{ \begin{array}{l} V_j(\Delta t)h_j^{n+1-} = V(0)h_j^n, \\ V_j(\Delta t)(hu)_j^{n+1-} = V(0)(hu)_j^n - V(0)\frac{\Delta t}{\Delta x}(p_{j+1/2}^* - p_{j-1/2}^*), \end{array} \right. $$

123

with $V(0) = \Delta x$ and according to $dV(t)/dt = V(t)\partial_x u(x(t; \xi), t)$,

124

$$ (14) \quad V_j(\Delta t) = V(0) + \int_0^{\Delta t} V(t) \partial_x u(x(t; \xi), t) dt \approx V(0) + \Delta t (u_{j+1/2}^* - u_{j-1/2}^*) = \Delta x L_j^n. $$

125

Therefore, one clearly sees that the discretizations (13) and (14) of (9) (or equivalently (11)) allow to recover (8). At last, note now that from (12),

$$ \int_0^{\Delta t} \partial_t x_{\pm}(t; x_{j\pm 1/2}) dt = \int_0^{\Delta t} u(x_{\pm}(t; x_{j\pm 1/2}), t) dt $$

so that

$$ x_{\pm}(\Delta t; x_{j\pm 1/2}) - x_{\pm}(0; x_{j\pm 1/2}) \approx \Delta t u_{j\pm 1/2}^*, $$

and by subtraction we also have

$$ V_j(\Delta t) := x_+(\Delta t; x_{j+1/2}) - x_-(\Delta t; x_{j-1/2}) = V(0) + \Delta t(u_{j+1/2}^* - u_{j-1/2}^*) = \Delta x L_j^n, $$

and therefore

$$ L_j^n = 1 + \frac{\Delta t}{\Delta x} (u_{j+1/2}^* - u_{j-1/2}^*) = \frac{x_+(\Delta t; x_{j+1/2}) - x_-(\Delta t; x_{j-1/2})}{\Delta x}. $$

The transport step. Let us now turn to the transport or remap step. In the framework of the present equivalent reformulation of our algorithm, we have now at hand some constant approximations denoted $h_j^{n+1-}$ and $(hu)_j^{n+1-}$ of $h$ and $(hu)$ at time $t^{n+1}$ on the new cells $(x_-(\Delta t; x_{j-1/2}), x_+(\Delta t; x_{j+1/2}))$ and whose volume is $V_j(\Delta t)$. In this context, the aim of the remap step is just to project these values on the original mesh, which simply writes with clear notations and $X = h, hu$,

$$ X_j^{n+1} = \frac{1}{\Delta x} \int_{x_{j-1/2}}^{x_{j+1/2}} X_j^{n+1-}(x) dx, $$
or equivalently

$$X_j^{n+1} = \frac{x_-(\Delta t; x_{j-1/2}) - x_{j-1/2}}{\Delta x} X_{j-1/2}^{n+1-} + \frac{x_+(\Delta t; x_{j+1/2}) - x_-(\Delta t; x_{j-1/2})}{\Delta x} X_j^{n+1-} + \frac{x_{j+1/2} - x_+(\Delta t; x_{j+1/2})}{\Delta x} X_{j+1/2}^{n+1-}$$

where we have set for all $j$

$$X_{j+1/2}^{n+1-} = \begin{cases} X_j^{n+1-} & \text{if } u_{j+1/2}^* \ge 0, \\ X_{j+1}^{n+1-} & \text{if } u_{j+1/2}^* < 0. \end{cases}$$

It is easy to check that these final values coincide with the ones of the previous formulation of the scheme.

**3. The Lagrange-Projection scheme with a non flat topography.** Now we would like to extend this approach to the case of a non flat bottom paying a particular attention to the well-balanced property. Here the Lagrange-Projection approach aims at solving successively the Lagrangian form of our system, namely

$$\left\{ \begin{aligned} \partial_t h + h \partial_x u &= 0, \\ \partial_t (hu) + h u \partial_x u + \partial_x \left( g \frac{h^2}{2} \right) + g h \partial_x z &= 0, \end{aligned} \right.$$

which gives in Lagrangian coordinates, with $\tau = 1/h$,

$$\left\{ \begin{aligned} \partial_t \tau - \partial_m u &= 0, \\ \partial_t u + \partial_m \left( g \frac{h^2}{2} \right) + g h \partial_m z &= 0, \end{aligned} \right.$$

and the transport system

$$\left\{ \begin{aligned} \partial_t h + u \partial_x h &= 0, \\ \partial_t (hu) + u \partial_x (hu) &= 0. \end{aligned} \right.$$

At this stage, it is important to have in mind that in order to get the fully well-balanced property, we are going to consider in-cell reconstructions in both Lagrangian and Transport steps. Moreover, these reconstructions will make both steps of the strategy to be fully well-balanced for $(hu)$, in the sense that if the solution at time $t^n$ is an equilibrium state, we will have $(hu)_j^{n+1} = (hu)_j^{n+1-} = (hu)_j^n$ for all $j$. This is clearly sufficient but not necessary as the fully well-balanced property simply asks for $(hu)_j^{n+1} = (hu)_j^n$ for all $j$. As far as $h$ is concerned, we will have $h_j^{n+1-} \neq h_j^n$ in general, but the overall algorithm will be fully well-balanced, namely $h_j^{n+1} = h_j^n$ for all $j$. Let us now describe the two steps in details.

**3.1. The Lagrangian step.** Considering first the Lagrangian step, following the ideas used in the flat topography case, we propose the following generalization

$$\left(15) \quad \begin{cases} L_j^n h_j^{n+1-} = h_j^n, \\ L_j^n (hu)_j^{n+1-} = (hu)_j^n - \frac{\Delta t}{\Delta x} (p_{j+1/2}^* - p_{j-1/2}^*) - \Delta t \{gh \partial_x z\}_j, \end{cases}$$

with

$$L_j^n = 1 + \frac{\Delta t}{\Delta x} (u_{j+1/2}^* - u_{j-1/2}^*),$$
and where

$$ \{gh\partial_x z\}_j \approx \frac{1}{\Delta x} \int_{x_{j-1/2}}^{x_{j+1/2}} gh(x, t^n)\partial_x z(x, t^n)dx. $$

At this stage, it remains to define $u_{j+1/2}^*$, $p_{j+1/2}^*$ and $\{gh\partial_x z\}_j$ for all $j$ in such a way that the fully well-balanced property holds true for this step.

**Reconstruction procedure.** In order to get this property, we propose to perform in-cell reconstructions of the variables at time $t^n$, that is to say just before starting the Lagrangian step. With this in mind, we suggest to understand the cell values $h_j^n$ and $(hu)_j^n$ as the averages of stationary solutions defined in the corresponding cell $C_j$. Thus, we aim at defining the functions $x \to h(x; K_{1,j}, K_{2,j})$ and $x \to u(x; K_{1,j}, K_{2,j})$ such that

$$ \left\{ \begin{array}{l} (hu)(x; K_{1,j}, K_{2,j}) = K_{1,j}, \\ \frac{u(x; K_{1,j}, K_{2,j})^2}{2} + g(h(x; K_{1,j}, K_{2,j}) + z(x)) = K_{2,j}, \end{array} \right. $$

with two constants $K_{1,j}$ and $K_{2,j}$ such that

$$ h_j^n = \frac{1}{\Delta x} \int_{x_{j-1/2}}^{x_{j+1/2}} h(x; K_{1,j}, K_{2,j})dx $$

and

$$ (hu)_j^n = \frac{1}{\Delta x} \int_{x_{j-1/2}}^{x_{j+1/2}} (hu)(x; K_{1,j}, K_{2,j})dx = K_{1,j}. $$

In practice, and since we are only interested in a first order scheme, it is more convenient to approximate the above integrals by the mid-point quadrature formula, which allows to identify the cell averages $h_j^n$ (respectively $(hu)_j^n$) and the mid-point values $h(x_j; K_{1,j}, K_{2,j})$ (resp. $x \to u(x_j; K_{1,j}, K_{2,j})$). We are thus led to define the functions $x \to h(x; K_{1,j}, K_{2,j})$ and $x \to u(x; K_{1,j}, K_{2,j})$ by

$$ \tag{16} \begin{array}{l} (hu)(x; K_{1,j}, K_{2,j}) = (hu)_j^n, \\ \frac{u(x; K_{1,j}, K_{2,j})^2}{2} + g(h(x; K_{1,j}, K_{2,j}) + z(x)) = \frac{(u_j^n)^2}{2} + g(h_j^n + z_j), \\ z_j = \frac{1}{\Delta x} \int_{x_{j-1/2}}^{x_{j+1/2}} z(x)dx \end{array} $$

for all $j$. Using the $\mathbf{U}$ notation, we have thus been able to reconstruct smooth stationary solutions $x \to \mathbf{U}(x; K_{1,j}, K_{2,j})$ of the SWE within each cell and the average of which coincides with the solution at time $t^n$ for the conservative variables and up to second-order accuracy. These stationary solutions are then used to define left and right extrapolated values at the mesh interfaces by setting for all $j$

$$ \mathbf{U}_{j+1/2-} = \mathbf{U}(x_{j+1/2}; K_{1,j}, K_{2,j}), \quad \mathbf{U}_{j-1/2+} = \mathbf{U}(x_{j-1/2}; K_{1,j}, K_{2,j}). $$

Remark that, to close the definition of the reconstructed values, the bottom topography should be defined at the interface. Following the ideas in [1] and [11], let us consider a reconstructed bottom at $z_{j+1/2}$, consistent with $z(x_{j+1/2})$, whose precise definition will be discussed later. Thanks to (16) and using the notation $q = hu$, we consider the reconstructed states $h_{j+1/2\pm}$, $q_{j+1/2\pm}$ given by

147

148

149

150

151
$$152 \quad (17) \quad \left\{ \begin{array}{l} q_{j+1/2+}^n = q_{j+1}^n, \\ \frac{(q_{j+1}^n)^2}{2(h_{j+1/2+}^n)^2} + g(h_{j+1/2+}^n + z_{j+1/2}) = \frac{(q_{j+1}^n)^2}{2(h_{j+1}^n)^2} + g(h_{j+1}^n + z_{j+1}), \end{array} \right.$$

and

$$153 \quad (18) \quad \left\{ \begin{array}{l} q_{j+1/2-}^n = q_j^n, \\ \frac{(q_j^n)^2}{2(h_{j+1/2-}^n)^2} + g(h_{j+1/2-}^n + z_{j+1/2}) = \frac{(q_j^n)^2}{2(h_j^n)^2} + g(h_j^n + z_j). \end{array} \right.$$

It is important to notice that the values $h_{i+1/2\pm}$ may not be uniquely defined or not defined at all. Let us recall that this idea of using in-cell reconstruction using equilibria has been used, for instance, in [28, 10, 11]. In particular, let us consider a similar notation and the following proposition in [10]:

PROPOSITION 3.1. Let $u_0 \in \mathbb{R}$, $h_0 \ge 0$, $\delta \in \mathbb{R}$ and define $q_0 = h_0 u_0$. Consider the system

$$161 \quad (19) \quad \left\{ \begin{array}{l} \frac{q_0^2}{2(h^*)^2} + gh^* = \frac{q_0^2}{2h_0^2} + g(h_0 + \delta), \\ h^*u^* = q_0, \\ h^* \ge 0, \ u^* \in \mathbb{R}, \end{array} \right.$$

and denote $h_s \equiv h_s(q_0) = \left(\frac{q_0^2}{g}\right)^{\frac{1}{3}}$. There exists a solution $(h^*, u^*)$ to (19) if and only if

$$162 \quad (20) \qquad \delta + \delta_0(h_0, q_0) \ge 0,$$

where

$$166 \qquad \delta_0(h_0, q_0) = g^{-1} \left( \frac{q_0^2}{2h_0^2} + gh_0 - \frac{3}{2} (g|q_0|)^{\frac{2}{3}} \right) \in \mathbb{R}^+.$$

Moreover,

i) If we have equality in (20), there is only one solution $(h^*, u^*)$ to (19), it is given by

$$169 \quad (21) \qquad h^* = h_s, \qquad u^* = \begin{cases} \frac{q_0}{h_s} & \text{if } q_0 \ne 0, \\ 0 & \text{if } q_0 = 0. \end{cases}$$

ii) If we have a strict inequality in (20), then there are exactly two different solutions $(h_{sup}^*, u_{sup}^*)$ and $(h_{sub}^*, u_{sub}^*)$ to (19), with $h_{sup}^* \le h_s < h_{sub}^*$, and $h_{sup}^* < h_s$ for $q_0 \ne 0$.

iii) A solution $(h^*, u^*)$ to (19), with $h^*u^* \ne 0$ satisfies:

• It is a critical point $((u^*)^2 = gh^*)$ if, and only if, $h^* = h_s$,

• It is a subcritical point $((u^*)^2 < gh^*)$ if, and only if, $h^* > h_s$,

• It is a supercritical point $((u^*)^2 > gh^*)$ if, and only if, $h^* < h_s$.

iv) Assume that $\delta$ is such that (20) is satisfied strictly so that there are exactly two solutions $h_{sub}$ and $h_{sup}$ for (19).
179

a) When $(h_0, u_0)$ is subcritical, we have

180

a.1) if $\delta > 0$, then $h_{sup}^* < h_0 < h_{sub}^*$,

181

a.2) if $\delta \le 0$, then $h_{sup}^* < h_{sub}^* \le h_0$,

182

b) When $(h_0, u_0)$ is supercritical, we have

183

b.1) if $\delta \ge 0$, then $h_{sup}^* \le h_0 < h_{sub}^*$,

184

b.2) if $\delta < 0$, then $h_0 \le h_{sup}^* < h_{sub}^*$,

185

c) When $(h_0, u_0)$ is critical, we have $h_{sup}^* \le h_0 < h_{sub}^*$.

186

□

187

Remark 3.2. Based on this result, solving equations (17) and (18) means:

188

a) When $\Delta z = z_{j+1} - z_j \to 0$, we can assure the existence of the reconstructed states. This means that when $U_j$ and $U_{j+1}$ are not critical points ($\delta_0(h_j, q_j) > 0$ and $\delta_0(h_{j+1}, q_{j+1}) > 0$) and the bottom is continuous, we can define the reconstructed values provided that the spatial resolution is high enough. Nevertheless, some problems may arise when we are near a critical point at the discrete level.

189

190

191

192

193

b) When the reconstructed states exist, we have in general two possible choices: a subcritical and supercritical value. We will select in general the solution that preserves the same character, that is, for $U_j$ subcritical (resp. supercritical) we choose $U_{j+1/2-}$ and $U_{j-1/2+}$ subcritical (resp. supercritical). This is true far from sonic points and for $\Delta z$ sufficiently small.

194

195

196

197

c) In the general case, the choice of $z_{j+1/2}$ should be taken with care. We require that

198

199

$$z_j - z_{j+1/2} \ge -\delta_0(h_j, q_j) \text{ and } z_{j+1} - z_{j+1/2} \ge -\delta_0(h_{j+1}, q_{j+1})$$

200

so that the reconstructed states exist. Moreover, we shall impose that the reconstructed states verify

201

202

$$h_{j+1/2-} \le h_j \text{ or } h_{j+1/2-} \le h_s(q_j),$$

203

and

204

$$h_{j+1/2+} \le h_{j+1} \text{ or } h_{j+1/2+} \le h_s(q_{j+1}),$$

205

which will play an essential role for the positivity of the scheme.

206

Based on those conditions we propose the following definition:

207

i) Case $z_{j+1} \ge z_j$:

208

• Case $z_{j+1} \le z_j + \delta_0(h_j, q_j)$:

209

We set $z_{j+1/2} = z_{j+1}$, $h_{j+1/2+} = h_{j+1}$ and $h_{j+1/2-}$ the solution of (18) with the same character (subcritical or supercritical) as $h_j$.

210

211

• Case $z_{j+1} > z_j + \delta_0(h_j, q_j)$:

212

We set $z_{j+1/2} = z_j + \delta_0(h_j, q_j)$, $h_{j+1/2-} = h_s(h_j, q_j)$ and $h_{j+1/2+}$ the supercritical solution of (17).

213

214

ii) Case $z_{j+1} < z_j$:

215

• Case $z_j \le z_{j+1} + \delta_0(h_{j+1}, q_{j+1})$:

216

We set $z_{j+1/2} = z_j$, $h_{j+1/2-} = h_j$ and $h_{j+1/2+}$ the solution of (17) with the same character (subcritical or supercritical) as $h_{j+1}$.

217

218

• Case $z_j > z_{j+1} + \delta_0(h_{j+1}, q_{j+1})$:

219

We set $z_{j+1/2} = z_{j+1} + \delta_0(h_{j+1}, q_{j+1})$, $h_{j+1/2-} = h_s(h_{j+1}, q_{j+1})$ and $h_{j+1/2+}$ the supercritical solution of (18).

220

d) Although the previous definitions have been written for a general case, remark that in most situations it becomes simpler: consider continuous bottom and assume a
spatial discretization sufficiently small, such that

$$z_{j+1} - z_j < \delta_0(h_j, q_j), \text{ for } z_{j+1} \ge z_j$$

and

$$z_j - z_{j+1} < \delta_0(h_{j+1}, q_{j+1}), \text{ for } z_j > z_{j+1}.$$

Then, we get $z_{j+1/2} = \max(z_j, z_{j+1})$ and the values $h_{j+1/2+}$ and $h_{j+1/2-}$ correspond to the solutions of (17) and (18) with the same character (subcritical or supercritical) as $h_{j+1}$ and $h_j$ respectively.

To conclude this paragraph, notice that the proposed reconstruction procedure can be understood as a generalization of the well-known hydrostatic reconstruction associated with the lake at rest equilibrium and defined by

$$ (22) \qquad \begin{cases} u(x; K_{1,j}, K_{2,j}) = u_j^n, \\ h(x; K_{1,j}, K_{2,j}) + z(x) = h_j^n + z_j, \end{cases} $$

see [1] for more details.

Now that the reconstruction is defined, we use the left and right traces $U_{j+1/2\pm}$ to calculate the numerical fluxes $u_{j+1/2}^*$ and $p_{j+1/2}^*$ by simply setting

$$u_{j+1/2}^* = u^*(U_{j+1/2-}, U_{j+1/2+}) \quad \text{and} \quad p_{j+1/2}^* = p^*(U_{j+1/2-}, U_{j+1/2+}),$$

and with

$$ (23) \qquad a_{i+1/2} = \max \left\{ h_{i+1/2+}^n \sqrt{gh_{i+1/2+}^n}, h_{i+1/2-}^n \sqrt{h_{i+1/2-}^n} \right\}. $$

Let us now define $\{gh\partial_x z\}_j$ such that

$$ \{gh\partial_x z\}_j \approx \frac{1}{\Delta x} \int_{x_{j-1/2}}^{x_{j+1/2}} gh(x, t^n) \partial_x z(x) dx. $$

Since we have reconstructed stationary solutions in each cell, we clearly have for all $j$

$$ \partial_x(hu^2 + g\frac{h^2}{2})(x; K_{1,j}, K_{2,j}) = -(gh\partial_x z)(x; K_{1,j}, K_{2,j}), $$

which, by integrating in the interior of the cell, gives

$$ \begin{aligned} -\{gh\partial_x z\}_j = & \frac{1}{\Delta x} (p(x_{j+1/2}; K_{1,j}, K_{2,j}) - p(x_{j-1/2}; K_{1,j}, K_{2,j})) + \\ & \frac{1}{\Delta x} (hu^2(x_{j+1/2}; K_{1,j}, K_{2,j}) - hu^2(x_{j-1/2}; K_{1,j}, K_{2,j})), \end{aligned} $$

that is to say

$$ -\{gh\partial_x z\}_j = \frac{p_{j+1/2-}^n - p_{j-1/2+}^n}{\Delta x} + \frac{(hu^2)_{j+1/2-}^n - (hu^2)_{j-1/2+}^n}{\Delta x}, $$

where of course $p_{j+1/2-}^n = g(h_{j+1/2-}^n)^2/2$ and $p_{j-1/2+}^n = g(h_{j-1/2+}^n)^2/2$. Since, $(hu)_{j+1/2-}^n = (hu)_{j-1/2+}^n = (hu)_{j}^n$, we finally get

$$ -\{gh\partial_x z\}_j = \frac{p_{j+1/2-}^n - p_{j-1/2+}^n}{\Delta x} + (hu)_j^n \frac{u_{j+1/2-}^n - u_{j-1/2+}^n}{\Delta x}. $$
230

Remark that the definition of {$gh∂_xz$}$_j$ is related to the integral of the source term along the path that corresponds to the integral curves defining stationary states. We refer the reader to [29] and [13] for further details.

231

232

233

**3.2. The Transport step.** Considering now the Transport step, we propose first to rewrite equation

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

in following way

$$\partial_t X + \partial_x(Xu) - X \partial_x u = 0,$$

with $X = h, (hu)$. Now we suggest the following update formulas

$$X_j^{n+1} = X_j^{n+1-} - \frac{\Delta t}{\Delta x} (X_{j+1/2}^{n+1-} u_{j+1/2}^* - X_{j-1/2}^{n+1-} u_{j-1/2}^*) + X_j^{n+1-} \frac{\Delta t}{\Delta x} (u_{j+1/2}^* - u_{j-1/2}^*),$$

234

or equivalently

$$(24) \qquad X_j^{n+1} = L_j^n X_j^{n+1-} - \frac{\Delta t}{\Delta x} (X_{j+1/2}^{n+1-} u_{j+1/2}^* - X_{j-1/2}^{n+1-} u_{j-1/2}^*),$$

235

236

where the $u_{j+1/2}^*$ have the same definition as in the Lagrangian step, and where the values $X_{j+1/2}^{n+1-}$ denote relevant approximations of $X = h, (hu)$ at the interface $x_{j+1/2}$ at time $t^{n+1-}$. Let us notice first that the whole scheme is actually conservative since easy calculations show that it writes

237

238

239

$$(25) \qquad \left\{ \begin{aligned} h_j^{n+1} &= h_j^n - \frac{\Delta t}{\Delta x} (h_{j+1/2}^{n+1-} u_{j+1/2}^* - h_{j-1/2}^{n+1-} u_{j-1/2}^*), \\ (hu)_j^{n+1} &= (hu)_j^n - \frac{\Delta t}{\Delta x} ((hu)_{j+1/2}^{n+1-} u_{j+1/2}^* + p_{j+1/2}^* - (hu)_{j-1/2}^{n+1-} u_{j-1/2}^* - \Delta t \{gh \partial_x z\}_j). \end{aligned} \right.$$

240

241

Finally, the values $X_{j+1/2}^{n+1-}$ have to be defined for $X = h, (hu)$ in order to get the expected well-balanced property. Let us now address this issue.

242

243

If one first considers the definition of $h_{j+1/2}^{n+1-}$, it is a matter of defining an approximation of $h$ at the interface $x_{j+1/2}$ and at time $t^{n+1-}$. Mimicking the definition used in the case of a flat topography, namely

$$h_{j+1/2}^{n+1-} = \begin{cases} h_j^{n+1-} & \text{if } u_{j+1/2}^* \ge 0, \\ h_{j+1}^{n+1-} & \text{if } u_{j+1/2}^* \le 0, \end{cases}$$

and the formula $L_j^n h_j^{n+1-} = h_j^n$ which was used to define an approximation of $h$ on the j-th cell and at time $t^{n+1-}$, we are tempted to set

$$h_{j+1/2}^{n+1-} = \begin{cases} h_{j+1/2-}^{n+1-} & \text{if } u_{j+1/2}^* \ge 0, \\ h_{j+1/2+}^{n+1-} & \text{if } u_{j+1/2}^* \le 0, \end{cases}$$

244

where the values $h_{j+1/2\pm}^{n+1-}$ have to be defined.

245

246

Definition of $h_{j+1/2\pm}^{n+1-}$. Here we propose to define

247

$$(26) \qquad h_{j+1/2\pm}^{n+1-} = h_{j+1/2\pm}^n$$
248 which is a simple way to account for the reconstruction procedure introduced in the
249 Lagrangian step. It will also guarantee the fully well-balanced property.

250 Note that following the same idea introduced in Section 2.3, we may interpret again
251 the transport step as the projection of the values in Lagrangian coordinates into the
252 Eulerian coordinates. The reconstructed states $h_{j+1/2\pm}^{n+1-}$ can thus be seen as the values
253 at time $t^{n+1}$ to the left and to the right of $x_{j+1/2}^* = x(\Delta t; x_{j+1/2})$, the position at time
254 $t^{n+1}$ of the moving Lagrangian intercell. Here we propose $h_{j+1/2\pm}^n$ as an approximation
255 of those values.

256
257 Regarding now the definition of $(hu)_{j+1/2}^{n+1-}$, it is a matter of defining again a suitable approximation of $(hu)$ at time $t^{n+1-}$. Here, we propose to keep the same definition used in the case of a flat topography, namely

$$
(hu)_{j+1/2}^{n+1-} = \begin{cases} (hu)_j^{n+1-} & \text{if } u_{j+1/2}^* \ge 0, \\ (hu)_{j+1}^{n+1-} & \text{if } u_{j+1/2}^* \le 0, \end{cases}
$$

258 which will be sufficient to get the fully well-balanced property, as we show now.

**3.3. Positivity and fully well-balanced property.** In this section, we prove that the Lagrangian-Projection scheme defined by (15)-(24) is fully well-balanced and that it preserves positivity of water height.

Well-balanced property. We shall check that the scheme defined previously preserves all smooth stationary states. Let us assume that the solution at time t^n corresponds to a general steady state of the SWE, that is to say

$$
(hu)_j^n = K_1 \quad \text{and} \quad \frac{(u_j^n)^2}{2} + g(h_j^n + z_j) = K_2, \quad \forall j \in \mathbb{Z},
$$

where $K_1$ and $K_2$ are two given constants. The reconstruction procedure gives $\mathbf{U}_{j+1/2-} = \mathbf{I}$
$\mathbf{U}_{j+1/2+}$ and in particular

$$
u_{j+1/2}^{*} = u_{j+1/2 \pm}^{n} := u_{j+1/2}^{n} \quad \text{and} \quad p_{j+1/2}^{*} = p_{j+1/2 \pm}^{n} := p_{j+1/2}^{n}.
$$

As a consequence, the Lagrangian step (15) gives

$$
\left\{
\begin{array}{l}
L_j^n h_j^{n+1-} = h_j^n, \\
L_j^n (hu)_j^{n+1-} = (hu)_j^n - \frac{\Delta t}{\Delta x} (p_{j+1/2}^n - p_{j-1/2}^n) - \Delta t \{gh \partial_x z\}_j,
\end{array}
\right.
$$

with

$$
\begin{align*}
L_j^n &= 1 + \frac{\Delta t}{\Delta x} (u_{j+1/2}^n - u_{j-1/2}^n), \\
-g \partial_x z \Big|_j &= \frac{p_{j+1/2}^n - p_{j-1/2}^n}{\Delta x} + (hu)_j^n \frac{u_{j+1/2}^n - u_{j-1/2}^n}{\Delta x},
\end{align*}
$$

that is to say

$$
\left\{
\begin{array}{l}
L_j^n h_j^{n+1-} = h_j^n, \\
L_j^n (hu)_j^{n+1-} = L_j^n (hu)_j^n.
\end{array}
\right.
$$

Under the usual CFL restriction, we thus have $(hu)_j^{n+1-} = (hu)_j^n$ so that the Lagrangian step is fully well-balanced for $(hu)$. Note that in general, $h_j^{n+1-} \neq h_j^n$, except if $u_{j+1/2} = 0$.
$u_{j-1/2}$ gets true, which happens in the particular case of the lake at rest steady state ($u_j^n = 0$ for all $j$).

We also have

$$h_{j+1/2-}^{n+1-} = h_{j+1/2-}^{n} = h_{j+1/2+}^{n} = h_{j+1/2+}^{n+1-}.$$

As a consequence, (25) gives

$$\left\{ \begin{array}{l} h_j^{n+1} = h_j^n - \frac{\Delta t}{\Delta x}((hu)_j^n - (hu)_j^n), \\ (hu)_j^{n+1} = (hu)_j^n - \frac{\Delta t}{\Delta x}((hu)_j^n u_{j+1/2}^n + p_{j+1/2}^n - (hu)_j^n u_{j-1/2}^n - p_{j-1/2}^n) \\ \qquad - \Delta t \{gh \partial_x z\}_j, \end{array} \right.$$

Therefore, the scheme is thus fully well-balanced since $h_j^{n+1} = h_j^n$ and $(hu)_j^{n+1} = (hu)_j^n$ for all $j$.

**Positivity property.** It is clear that the Lagrangian step (15) preserves positivity of water thickness provided that $L_j^n > 0$, which will be true under a suitable CFL condition. The question is then whether the transport step preserves positivity or not. Following the ideas presented in [1], given a numerical scheme for water height in the form

$$h_j^{n+1} = h_j^n - \frac{\Delta t}{\Delta x} (\mathcal{F}_{j+1/2}^h - \mathcal{F}_{j-1/2}^h),$$
the scheme will preserve the positivity of $h$ if whenever $h_j^n = 0$, one has

$$\mathcal{F}_{j+1/2}^h - \mathcal{F}_{j-1/2}^h \le 0.$$

The fully well-balanced Lagrange-Projection for variable $h$ writes

$$h_j^{n+1} = h_j^{n+1-} - \frac{\Delta t}{\Delta x} \left( u_{j+1/2}^* (h_{j+1/2}^{n+1-} - h_j^{n+1-}) + u_{j-1/2}^* (h_j^{n+1-} - h_{j-1/2}^{n+1-}) \right).$$

Now, definitions given in Remark 3.2. c) grants that when $U_j^n$ is a subcritical state, then

$$\max(h_{j+1/2-}^{n+1-}, h_{j-1/2+}^{n+1-}) \le h_j^n.$$

Conversely, when $U_j^n$ is supercritical, we have

$$\max(h_{j+1/2-}^{n+1-}, h_{j-1/2+}^{n+1-}) \le h_s(q_j^n).$$

It follows then that whenever $h_j^n = 0$ we get

$$h_j^{n+1-} = 0, \quad h_{j+1/2-}^{n+1-} = 0, \quad \text{and} \quad h_{j-1/2-}^{n+1-} = 0.$$

In that case, we get

$$h_j^{n+1} = -\frac{\Delta t}{\Delta x} \left( u_{j+1/2}^* h_{j+1/2}^{n+1-} - u_{j-1/2}^* h_{j-1/2}^{n+1-} \right).$$

Now, it is easy to check that in this case

$$u_{j+1/2}^* h_{j+1/2}^{n+1-} = (u_{j+1/2}^*)^{-} h_{j+1/2+}^{n+1-} \le 0,$$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Subcritical steady state</figcaption></figure>

and

$$u_{j-1/2}^{*}h_{j-1/2}^{n+1-} = (u_{j-1/2}^{*})^{+}h_{j-1/2+}^{n+1-} \geq 0,$$

and the result follows.

## 4. Numerical results.

### 4.1. Subcritical steady state.
In this first test we will verify the fully well balanced property of the scheme. To do so, let us consider a bottom topography

$$ (27) \qquad z(x) = 0.5 \exp \left( -\frac{(x - 0.5)^2}{0.005} \right), $$

and the subcritical steady state solution verifying

$$ hu = 0.5, \quad \frac{u^2}{2} + g(h + z) = 0.125. $$

The simulation is carried out with 200 cells and open boundary conditions. As expected, this steady state is preserved by the scheme with an error $10^{-12}$. Figure 1 shows the surface elevation as well as the discharge at time $t = 1$ which corresponds just to the given initial condition.

### 4.2. Generation of subcritical steady state.
The objective of this test is to study the convergence of the scheme to steady states. To do so let us consider the bottom topography (27) and we consider as initial condition a lake at rest situation: $h + z = 1$, $u = 0$. We impose at the left boundary the condition $hu(x = 0, t) = 0.5$ and at the right boundary the condition $h(x = 1, t) = 1$. It is expected that, for long time, the solution will converge to the steady state presented previously. The simulation is carried out with 200 cells and the time evolution of the surface and discharge are shown in Figure 2. Here a comparison with the solutions computed with 800 cells are also shown. Figure 3 shows the steady state reached for $t = 50$ and comparison with the exact reference solution. Moreover, $10^{-12}$. Figure 4 shows the deviation from the expected constant values for $q$ and $\frac{u^2}{2} + g(h + z)$ at the steady state. We see that the
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Generation of subcritical steady state from stationary initial condition: surface and discharge evolution. Continuous lines correspond to the solution computed with 200 cells and discontinuous to the one computed with 800 cells</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Generation of subcritical steady state from stationary initial condition: surface elevation compared to reference solution</figcaption></figure>

achieved error is similar to the previous case, order $10^{-12}$ for the first one and $10^{-12}$ for the second.

**4.3. Transcritical continuous steady state.** As it is known, the generalized hydrostatic reconstruction procedure may present some problems at critical points where the reconstructed states may not be well defined. We propose now to take as initial condition a transcritical steady state to check how this situation is handled by the scheme. We consider again the bottom topography (27) and we solve numerically the relations

$$hu = 1.5, \quad \frac{u^2}{2} + g(h + z) = E,$$
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Generation of subcritical steady state from stationary initial condition: deviation from the expected constant state</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Transcritical steady state: surface elevation</figcaption></figure>

with

$$E = \frac{1.5}{2h_s^2} + g(h_s + 0.5), \quad \text{and} \quad h_s = \left(\frac{1.5^2}{g}\right)^{1/3}.$$

We consider the subcritical solution for $x < 0.5$ and the supercritical solution for $x > 0.5$. Previous values guarantee that these solutions exist and there is a critical point at $x = 0.5$.

Figure 5 shows the surface elevation for time $t = 1$ using 201 cells. We remark that the scheme behaves well at $x = 0.5$ where the critical point is located. Let us recall that critical points may be problematic for the reconstruction using equilibria. This is due to the fact that at discrete level we represent the exact smooth functions by a piece-wise constant function consisting on the cell averages which gives in practice a jump on the bottom at the interfaces. Nevertheless we check in Figure 5(b) that the scheme captures correctly the solution. Figure 6 shows the deviation from the expected constant values for $q$ and $\frac{u^2}{2} + g(h+z)$ at the steady state. We see that the error is of order $10^{-14}$ for the first one and $10^{-13}$ for the second which shows that the steady state is preserved.
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Transcritical steady state: deviation from the expected constant state</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Perturbation of a steady state at times $t = 0, 0.02$ and $0.04$</figcaption></figure>

**4.4. Small perturbation of a steady state.** We consider now a small perturbation of previous transcritical steady state. To do so, let us consider the same bottom and initial condition and we add a perturbation in the water thickness by increasing its value in 0.15 for $0.15 < x < 0.25$.

Figures 7 and 8 show the evolution of the surface and the discharge at the beginning of the simulation, where the perturbation propagates along the domain. Figure 9 shows the surface elevation and discharge at time $t = 10$ when the steady state is reached. Compared to the prescribed initial steady state, we remark that the perturbation has not affected the reached steady state. Nevertheless, a small artifact is seen at the critical point in Figure 9(b) although it is very small. This small artifact reduces as the number of cells is increased.

Consider now a smaller perturbation of 0.0001 for $0.15 < x < 0.2$ and let us compare the scheme provided here with a similar based on just an hydrostatic reconstruction [1]. While the former is fully well-balanced, the latter is only well-balanced for lake at rest steady states. The simulations are done using 3201 volume cells. Figure 10 shows at time $t = 0.05$ the difference between the computed values $h$ and $q$ with respect to the unperturbed steady-state $h_0$ and $q_0$ described in Subsection 4.3. As we see, just from the beginning, the non fully well-balanced scheme is producing a perturbation in the
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Perturbation of a steady state at times t = 0.06, 0.1 and 0.12</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9. Perturbation of a steady state at times t = 10</figcaption></figure>

middle of the domain which is even bigger than the perturbation introduced on the initial condition. Nevertheless, the fully well-balanced scheme introduced here does not present this problem.

**4.5. Transcritical steady state with shock.** We consider now a numerical test taken from [9] consisting of a transcritical flow with a shock over a bump. The interval is [0, 25], the initial condition and depth bottom is given by

$$ (28) \qquad \begin{array}{ll} h(x, t = 0) = 0.33, & hu(x, t = 0) = 0.18, \\ z(x) = \begin{cases} 0.2 + 0.05(x - 10)^2, & \text{if } 8 < x < 12, \\ 0, & \text{otherwise,} \end{cases} & \end{array} $$

and the boundary conditions are given by $hu(x = 0, t) = 0.18$, $h(x = 25, t) = 0.33$.

We consider a uniform mesh composed by 201 cells.

The initial condition will evolve until a transcritical solution with a stationary shock is developed. Figure 11 shows the surface and discharge evolution obtained with the scheme proposed here. In Figure 12 we see the surface profile obtained and a comparison for a reference solution with increasing number of cells. We remark there is a small problem for the scheme at the critical point near $x = 10$ and the position of the shock
<figure><img src="image_11.png" /><figcaption>FIGURE 10. Perturbation of a steady state. Comparison with the unperturbed steady state at time $t = 0.05$ for fully and non fully well-balanced schemes</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 11. Transcritical steady state with shock: surface elevation and discharge for times $t=0,5,10,50$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 12. Transcritical steady state with shock: surface elevation</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIGURE 13. Transcritical steady state with shock: steady state condition for different number of cells</figcaption></figure>

<table><caption>TABLE 1
$L^1$ errors and numerical orders of accuracy.</caption><thead><tr><th rowspan="2">No. of cells</th><th colspan="2">h</th><th colspan="2">hu</th></tr><tr><th>L<sup>1</sup> error</th><th>order</th><th>L<sup>1</sup> error</th><th>order</th></tr></thead><tbody><tr><td>25</td><td>2.31e-02</td><td>0.00</td><td>1.53e-01</td><td>0.00</td></tr><tr><td>50</td><td>1.77e-02</td><td>0.38</td><td>1.09e-01</td><td>0.48</td></tr><tr><td>100</td><td>1.18e-02</td><td>0.59</td><td>6.67e-02</td><td>0.71</td></tr><tr><td>200</td><td>7.00e-03</td><td>0.75</td><td>3.62e-02</td><td>0.88</td></tr><tr><td>400</td><td>3.74e-03</td><td>0.91</td><td>1.79e-02</td><td>1.01</td></tr></tbody></table>

does not exactly coincide with reference solution. Nevertheless those problems disappear when the mesh is refined and the numerical solutions converge to the reference solution. Figure 13 shows the steady state conditions and comparison with a reference solution for different number of cells.

**4.6. Formal convergence of the scheme.** We intend now to study the formal convergence of the scheme proposed here. To do so, let us consider the initial condition given by

$$z(x) = 0.1 \cos(2\pi x), \quad h(x, 0) = 1.1 + 0.1 \sin(4\pi x) - z(x), \quad hu(x, 0) = 0.0.$$

We compute a reference solution using 3200 points in the interval $[0, 1]$ and using periodic boundary conditions up to time $t = 0.2$. Then we make a comparison with the solutions given by different number of cells. The $L^1$ errors are shown in Table 1. As expected, the scheme presented here is a first order scheme.

**4.7. Traveling wave over an obstacle.** We consider now bottom topography and initial condition given by

$$(29) \quad z(x) = \begin{cases} 0.5, & \text{if } 0.6 < x < 0.7, \\ 0, & \text{otherwise,} \end{cases} \quad h(x, t=0) = \begin{cases} 0.5, & \text{if } 0.6 < x < 0.7, \\ 1, & \text{otherwise,} \end{cases}$$

and $hu(x, t = 0) = 0$. This will produce a rarefaction followed by a shock that travels towards the obstacle at the bottom. We compute the solution using 201 cells in the
<figure><img src="image_15.png" /><figcaption>FIGURE 14. *Travelling wave over an obstacle: surface elevation and discharge for times t = 0, 0.04, 0.12, 0.2*</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 15. *Travelling wave over an obstacle: surface elevation and discharge for time t = 1*</figcaption></figure>

interval [0, 1] with open boundary conditions and compare the results with a reference solution computed on a fine grid with 3201 points. The time evolution of the surface and discharge are shown in Figure 14 where dashed lines correspond to the results using the finer grid and continuous line to the coarser grid. The final time is shown in 15 once the waves have left the domain and the steady state is reached. We see that scheme is able to reproduce the evolution of the waves in the flat region as well as when the shock encounters the obstacle.

**5. Conclusions.** We have introduced a positive and fully well-balanced Lagrange-Projection scheme for the Shallow Water Equations (SWE). This means that the scheme is able to preserve stationary solutions of the model with non zero velocity. Moreover, the scheme is positive under suitable CFL condition. Numerical results have shown that the scheme behaves well even in transcritical regimes or near critical points, provided that a sufficient number of cells is used. In a next work we propose to define semi-implicit fully well-balanced numerical schemes for SWE. Fully well-balanced high order Lagrange-Projection schemes will be also considered.

**Acknowledgements.** This research has been partially supported by the Spanish Government and FEDER through the coordinated Research projects MTM 2015-70490-
387

C2-1-R and MTM 2015-70490-C2-2-R. The authors were partially supported by the In-
novative Training Networks (ITN) grant 642768 (ModCompShock). The second author
was partially supported by the ANR project ACHYLLES (ANR-14-CE25-0001-03).

388

389

390

REFERENCES

391

[1] E. Audusse, F. Bouchut, M. O. Bristeau, R. Klein, B. Perthame, *A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows*, SIAM J. Sci. Comp., 25, 2050–2065, 2004.

392

393

394

[2] E. Audusse, C. Chalons, P. Ung, *A very simple well-balanced positive and entropy-satisfying scheme for the shallow-water equations*, to appear in Commun. in Math. Sci, 2015.

395

[3] A. Bermudez, M.E. Vazquez-Cendon, *Upwind Methods for Hyperbolic Conservation Laws with Source Terms*, Comp. & Fluids, 23, 1049–1071, 1994.

397

398

[4] C. Berthon, C. Chalons, *A fully well-balanced, positive and entropy-satisfying Godunov-type method for the Shallow-Water Equations*, Math. of Comp., 85, 1281–1307, 2016.

399

400

[5] C. Berthon, C. Chalons, S. Cornet, and G. Sperone, *Fully well-balanced, positive and simple approximate Riemann solver for shallow-water equations*, Bull. Braz. Math. Soc., New Series, 47(1), 1–14, 2016.

401

402

[6] C. Berthon, F. Foucher, *Efficient wellbalanced hydrostatic upwind schemes for shallow water equations*, J. Comput. Phys., 231, 4993–5015, 2012.

403

404

[7] G. Bispen, K. R. Arun, M. Lukacova-Medvidova, S. Noelle, *IMEX large time step finite volume methods for low Froude number shallow water flows*, Communications in Computational Physics, 16, 307–347, 2014.

405

406

[8] F. Bouchut, *A reduced stability condition for nonlinear relaxation to conservative laws*, J. Hyp. Diff. Eq., 1(1), 149–170, 2004.

407

408

[9] F. Bouchut, *Non-linear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources*, Frontiers in Mathematics, Birkhauser, 2004.

409

410

[10] F. Bouchut, T. Morales, *A subsonic-well-balanced reconstruction scheme for shallow water flows*, Siam J. Numer. Anal., 48(5), 1733–1758, 2010.

411

412

[11] M. J. Castro, A. Pardo, C. Parés, *Well-balanced numerical schemes based on a generalized hydrostatic reconstruction technique*, Mathematical Models and Methods in Applied Sciences, 17, 2065–2113, 2007.

413

414

[12] M.J. Castro, J.A. López, C. Parés. *High order exactly well-balanced numerical methods for shallow water systems*. Journal of Computational Physics, 246:242–264, 2013.

415

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

416

417

[14] C. Chalons, M. Girardin, and S. Kokh, *Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms*, SIAM J. Sci. Comput., 35(6), pp. a2874–a2902, 2013.

418

419

[15] C. Chalons, M. Girardin, and S. Kokh, *Operator-splitting based AP schemes for the 1D and 2D gas dynamics equations with stiff sources*, AIMS Series on Applied Mathematics, vol 8, 607–614, 2014.

420

421

422

[16] C. Chalons, M. Girardin, and S. Kokh, *An all-regime Lagrange-Projection like scheme for the gas dynamics equations on unstructured meshes*, Communications in Computational Physics (CICP), Vol. 20, No. 1, pp. 188-233, 2016.

423

424

[17] C. Chalons, P. Kestener, S. Kokh and M. Stauffert, *A large time-step and well-balanced Lagrange-Projection type scheme for the shallow-water equations*, to appear in Communications in Mathematical Sciences, 2017.

425

426

[18] G.-Q. Chen, C.D. Levermore, T.-P. Liu, *Hyperbolic conservation laws with stiff relaxation terms and entropy*, Comm. Pure Appl. Math., volume 47(6), 787–830, 1994.

427

428

429

[19] A. Chinnayya, A.-Y. LeRoux, N. Seguin, *A well-balanced numerical scheme for the approximation of the shallow-water equations with topography: the resonance phenomenon*, International Journal on Finite Volume (electronic), volume 1(1), 1–33, 2004.

430

431

[20] F. Coquel, Q.L. Nguyen, M. Postel, and Q.H. Tran, *Entropy-satisfying relaxation method with large time-steps for Euler IBVPs*, Math. Comp., 79, 1493–1533, 2010.

432

433

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

434

435

[22] G. Gallice, *Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source*, C. R. Math. Acad. Sci. Paris 334, no. 8, 713–716, 2002.

436

437

[23] G. Gallice, *Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equa-
446 tions in Lagrangian or Eulerian coordinates, Numer. Math. 94 no.4, 673-713, 2003.

447 [24] E. Godlewski and P.-A. Raviart, *Numerical approximation of hyperbolic systems of conservation laws*, Applied Mathematical Sciences, vol. 118, Springer-Verlag, New York, 1996.

448 [25] L. Gosse, *A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms*, Comput. Math. Appl., 39, 135-159, 2000.

449 [26] L. Gosse, *Computing qualitatively correct approximations of balance laws. Exponential-fit, well-balanced and asymptotic-preserving*, SEMA SIMAI Springer Series 2, 2013.

450 [27] J. M. Greenberg, A. Y. Leroux, *A well-balanced scheme for the numerical processing of source terms in hyperbolic equations*, SIAM J. Numer. Anal., 33, 1-16, 1996.

451 [28] S. Noelle, Y. Xing, and C.-W. Shu, *High-order well-balanced finite volume WENO schemes for shallow water equation with moving water*, Journal of Computational Physics, 226(1):29-58, 2007.

452 [29] Parés, C. *Numerical methods for nonconservative hyperbolic systems: a theoretical framework*, SIAM J. Numer. Anal., 44, 300-321, 2006.

453 [30] B. Perthame, C. Simeoni, *A kinetic scheme for the Saint-Venant system with source term*, Calcolo, 38, Springer-Verlag, 201-231, 2001.

454 [31] G. Russo, A. Khe. *High Order Well-Balanced Schemes Based on Numerical Reconstruction of the Equilibrium Variables*. Proceedings WASCOM 2009, Waves and Stability in Continuous Media, Palermo, June 28th to July 1st, 2009.

455 [32] I. Suliciu, *On the thermodynamics of fluids with relaxation and phase transitions*, Fluids with relaxation. Int. J. Engag. Sci., 36, 921-947, 1998.

456 [33] Y. Xing, *Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium*, J. Comput. Phys., 257, 536-553, 2014.

457 [34] Y. Xing, C.-W. Shu, *High order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms*, J. Comput. Phys., 214, 567-598, 2006.

458 [35] Y. Xing, C.-W. Shu, S. Noelle, *On the advantage of well-balanced schemes for moving-water equilibria of the shallow water equations*, J. Sci. Comput., 48, 339-349, 2011.

459 [36] Y. Xing, X. Zhang, C.-W. Shu, *Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations*, Adv. Water Resour., 33, 1476-1493, 2010.

460 [37] H. Zakerzadeh, *On the Mach-uniformity of the Lagrange-Projection scheme*, submitted, 2016.